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The density-matrix renormalization group (DMRG) is a numerical algorithm for the efficient trun- 
cation of the Hilbert space of low-dimensional strongly correlated quantum systems based on a 
rather general decimation prescription. This algorithm has achieved unprecedented precision in 
the description of one-dimensional quantum systems. It has therefore quickly acquired the status 
of method of choice for numerical studies of one-dimensional quantum systems. Its applications 
to the calculation of static, dynamic and thermodynamic quantities in such systems are reviewed. 
The potential of DMRG applications in the fields of two-dimensional quantum systems, quan- 
tum chemistry, three-dimensional small grains, nuclear physics, equilibrium and non-equilibrium 
statistical physics, and time-dependent phenomena is discussed. This review also considers the 
theoretical foundations of the method, examining its relationship to matrix-product states and 
the quantum information content of the density matrices generated by DMRG. 
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Consider a crystalline solid; it consists of some 10^^ 
or more atomic nuclei and electrons for objects on a hu- 
man scale. All nuclei and electrons are subject to the 
strong, long-range Coulomb interaction. While this sys- 
tem should a priori be described by the Schrodinger equa- 
tion, its explicit solution is impossible to find. Yet, it 
has become clear over the decades that in many cases 
the physical properties of solids can be understood to a 
very good approximation in the framework of some ef- 
fective one-body problem. This is a consequence of the 
very efficient interplay of nuclei and electrons to screen 
the Coulomb interaction but on the very shortest length 
scales. 

This comforting picture may break down as various ef- 
fects invalidate the fundamental tenet of weak effective 
interactions, taking us into the field of strongly corre- 
lated quantum systems, where the full electronic many- 
body problem has to be considered. Of the routes to 
strong correlation, let me mention just two, because it 
is at their meeting point that, from the point of view of 
physical phenomena, the density-matrix renormalization 
group (DMRG) has its origin. 

On the one hand, for arbitrary interactions, the Fermi 
liquid picture breaks down in one-dimensional solids and 
gives way to the Tomonaga-Luttinger liquid picture: es- 
sentially due to the small size of phase-space, scattering 
processes for particles close to the Fermi energy become 
so important that the (fermionic) quasiparticle picture 
of Fermi liquid theory is replaced by (bosonic) collective 
excitations. 

On the other hand, for arbitrary dimensions, screening 
is typically so effective because of strong delocalization 
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of valence electrons over the lattice. In transition metals 
and rare earths, however, the valence orbitals are inner 
d- or /-orbitals. Hence valence electrons are much more 
localized, albeit not immobile. There is now a strong 
energy penalty to place two electrons in the same local 
valence orbital, and the motion of valence electrons be- 
comes strongly correlated on the lattice. 

To study strongly correlated systems, simplified model 
Hamiltonians that try to retain just the core ingredients 
needed to describe some physical phenomenon and meth- 
ods for their treatment have been designed. Localization 
suggests the use of tight-binding lattice models, where lo- 
cal orbitals on one site can take iVgitc — 4 different states 
of up to two electrons (|0), 1 1), | i), | Ti>)- 

The simplest model Hamiltonian for just one va- 
lence orbital (band) with a kinetic energy term (elec- 
tron hopping between sites i and j with amplitude ty) 
and Coulomb repulsion is the on-site Hubbard model 
llfubbard, 1963, 1964a b), where just the leading on-site 
Coulomb repulsion U has been retained: 

-^Hubbard - ^ ti.i{cl^Cj^+h.C.)+U'^ni^n^l. (1) 

designates bonds. In the limit U/t ^ 1 double occu- 
pancy I ti) can be excluded resulting in the A^sitc = 3 
state t-J model: 

where the spin-spin interaction Jij — Atfj/U is due to 
a second-order virtual hopping process possible only for 
electrons of opposite spin on sites i and j. At half- filling, 
the model simplifies even further to the spin-i isotropic 
Heisenberg model, 

-^Heisenberg — ^ ^ Jij^i ' Sj" , (3) 
{ij) 

placing collective (anti)ferromagnetism, which it de- 
scribes, into the framework of strongly correlated sys- 
tems. These and other model Hamiltonians have been 
modified in multiple ways. 

In various fields of condensed matter physics, such 
as high-temperature superconductivity, low-dimensional 
magnetism (spin chains and spin ladders), low- 
dimensional conductors, and polymer physics, theoretical 
research is focussing on the highly nontrivial properties 
of these seemingly simple models, which are believed to 
capture some of the essential physics. Recent progress 
in ex periments on ultracold atomic gases l)Greiner et al\ . 
120021) has allowed the preparation of strongly correlated 
bosonic systems in optical lattices with tunable interac- 
tion parameters, attracting many solid-state physicists to 
this field. On the conceptual side, important old and new 
questions relating e.g. to topological effects in quantum 
mechanics, the wide field of quantum phase transitions, 
and the search for exotic forms of order stabilized by 



quantum effects at low temperatures are at the center of 
the physics of strongly correlated quantum systems. 

The typical absence of a dominant, exactly solvable 
contribution to the Hamiltonian, about which a per- 
turbative expansion such as in conventional many-body 
physics might be attempted, goes a long way in explain- 
ing the inherent complexity of strongly correlated sys- 
tems. This is why, apart from some exact solutions such 
as provided by the Bethe ansatz or certain toy models, 
most analytical approaches are quite uncontrolled in their 
reliability. While these approaches may yield important 
insight into the nature of the physical phenomena ob- 
served, it is ultimately up to numerical approaches to 
assess the validity of analytical approximations. 

Standard methods in the field are the exact diagonal- 
ization of small systems, the various quantum Monte 
Carlo methods, and, less widely used, coupled cluster 
methods and series expans ion techniques. Since its incep- 
tion by Steve lWhitd l)l992() . the density-matrix renormal- 
ization group has quickly achieved the status of a highly 
reliable, precise, and versatile numerical method in the 
field. Its main advantages are the ability to treat unusu- 
ally large systems at high precision at or very close to zero 
temperature and the absence of the negative sign problem 
that plagues the otherwise most powerful method, quan- 
tum Monte Carlo, making the latter of limited use for 
frustrated spin or fermionic systems. A first striking illus- 
tration of DMRG precision was given by White and Hu^ 
l|l993|) . Using modest numerical means to study the 
S = 1 isotropic antiferromagnetic Heisenberg chain, they 
calculated the ground state energy at almost machine 
precision, Eq = -1.401484038971(4) J, and the (Hal- 
dane) excitation gap as A = 0.41050(2) J. The spin-spin 
correlation length was found to be ^ = 6.03(2) lattice 
spacings. More examples of the versatility and precision 
which have made the reputation of the DMRG method 
will be found throughout the text. The major drawback 
of DMRG is that it displays its full force mainly for one- 
dimensional systems; nevertheless, interesting forays into 
higher dimensions have been made. By now, DMRG has, 
despite its complexity, become a standard tool of compu- 
tational physics that many research groups in condensed 
matter physics will want to have at hand. In the simula- 
tion of model Hamiltonians, it ranks to me as part of a 
methodological trias consisting of exact diagonalization, 
quantum Monte Carlo and DMRG. 

Most DMRG applications still treat low-dimensional 
strongly correlated model Hamiltonians, for which the 
method was originally developed. However, the core idea 
of DMRG, the construction of a renormalization-group 
fiow in a space of suitably chosen density matrices, is 
quite general. In fact, another excellent way of concep- 
tual thinking about DMRG is in the very general terms of 
quantum information theory. The versatility of the core 
idea has allowed the extension of DMRG applications to 
fields quite far away from its origins, such as the physics 
of (three-dimensional) small grains, equilibrium and far- 
from-equilibrium problems in classical and quantum sta- 



3 



tistical mechanics as well as the quantum chemistry of 
small to medium-sized molecules. Profound connections 
to exactly solvable models in statistical mechanics have 
emerged. 

The aim of this review is to provide both the algorith- 
mic foundations of DMRG and an overview of its main, 
by now quite diverse fields of applications. 

I start with a discussion of the key algorithmic ideas 
needed to deal with the most conventional DMRG prob- 
lem, the study of T = static properties of a one- 
dimensional quantum Hamiltonian (Sec. ^J- This in- 
cludes standard improvements to the basic algorithm 
that should always be used to obtain a competitive 
DMRG application, a discussion of the assessment of 
the numerical quality of DMRG output, as well as an 
overview of applications. Having read this first major 
section, the reader should be able to set up standard 
DMRG having the major algorithmic design problems in 
mind. 

I move on to a section on DMRG theory (Sec. 
discussing the properties of the quantum states DMRG 
generates (matrix product states) and the properties of 
the density matrices that are essential for its success; the 
section is closed by a reexamination of DMRG from a 
quantum information theory point of view. It might be 
quite useful to have at least some superficial grasp of the 
key results of this section in order to best appreciate the 
remainder of the review, while a more thorough reading 
might serve as a wrap-up in the end. All the sections 
that come after these key conceptual sections can then, 
to some large degree, be read independently. 

Among the various branches of applied DMRG, the 
applications to dynamical properties of quantum systems 
are presented first (Sec.lIV|l. Next, I discuss attempts to 
use DMRG for systems with potentially large numbers 
of local degrees of freedom such as phononic models or 
Bose-Hubbard models fSeclvTl. 

Moving beyond the world of T = physics in one di- 
mension, DMRG as applied to two-dimensional quantum 
Hamiltonians in real space with short-ranged interactions 
is introduced and various algorithmic variants are pre- 
sented (SecEB- 

Major progress has been made by abandoning the con- 
cept of an underlying real-space lattice, as various au- 
thors have developed DMRG variants for momentum- 
space DMRG (Sec. IVII.A^, DM RG for small molecule 
quantum chemistry fSec. IVII.BI) . and DMRG for meso- 
scopic small grains and nuclei (Sec. IVII.G|) . All these 
fields are currently under very active development. 

Early in the history of DMRG it was realized that its 
core renormalization idea might also be applied to the 
renormalization of transfer matrices that describe two- 
dimensional classical (Sec. IVIILAjl or one-dimensional 
quantum systems (Sec. IVIII.G|) in equilibrium at finite 
temperature. Here, algorithmic details undergo major 
changes, such that this class of DMRG methods is often 
referred to as TMRG (transfer matrix renormalization 
group ). 



Yet another step can be taken by moving to physical 
systems that are out of equihbrium fSec. llX|l . DMRG has 
been successfully applied to the steady states of nonequi- 
librium systems, leading to a non-hermitian extension of 
the DMRG where transition matrices replace Hamiltoni- 
ans. Various methods for time-dependent problems are 
under development and have recently started to yield re- 
sults unattainable by other numerical methods not only 
at T = 0, but also for finite temperature and in the pres- 
ence of dissipation. 

In this review, details of computer implementation 
have been excluded. Rather, I have tried to focus on 
details of algorithmic structure and their relation to the 
physical questions to be studied using DMRG, and to 
give the reader some idea about the power and the lim- 
itations of the method. In the more standard fields of 
DMRG, I have not been (able to be) exhaustive even in 
merely listing applications. In the less established fields 
of DMRG, I have tried to be more comprehensive and 
to provide as much discussion of their details as pos- 
sible, in order to make these fields better known and, 
hopefully, to stimulate thinking about further algorith- 
mic progress. I have excluded, for lack of space, any 
extensive discussions of the analytical methods whose 
development has been stimulated by DMRG concepts. 
Last but not least it should be mentioned that some of 
the topi cs of th i s revi ew have been considered by other 
authors: IWhitel l)l998() gives an introduction to the fun- 
damentals of DMRG; a very detailed survey of DMRG 
as it was in late 1998 has been provided in a collection of 
lectures and articles l|Peschel et a^J . IT999a|) : it also con- 
tains an account of DMRG history by White. More re- 
cently, the application of TMRG to quantum syst ems and 
two-d i mensional DMRG have been revieved bv IShibatal 
l|2on,^^ ■ lHallb'^^l2on.1^ gives a rather complete ov erview 
of DMRG applications. iDukelskv and PitteJ l|2004tl focus 
on DMRG applications to finite Fermi systems such as 
small grains, small molecules and nuclei. 

A word on notation: All state spaces considered here 
can be factorized into local state spaces {|(t)} labelled 
by Greek letters. DMRG forms blocks of lattice sites; 
(basis) states |m) of such blocks I denote by Latin let- 
ters. These states depend on the size of the block; when 
necessary, I indicate the block length by subscripts \mi). 
Correspondingly, \ai) is a local state on site i. Moreover, 
DMRG typically operates with two blocks and two sites, 
which are referred to as belonging to "system" or "envi- 
ronment" . Where this distinction matters, it is indicated 
by superscripts, \m^) or \a^). 



II. KEY ASPECTS OF DMRG 

Historically, DMRG has its origin in the analysis 
bv IWhite and Noack " ('1992') of the failure of real-space 
renormalization group (RSRG) methods to yield quanti- 
tatively acceptable results for the low-energy properties 
of quantum many-body problems. Most of the DMRG al- 
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gorithm can be formulated in standard (real space) renor- 
malization group language. Alternative points of view in 
terms of matrix product states and quantum information 
theory will be taken later fSec. lIII. Al and Sec. lIII.CI) . Be- 
fore moving on to the details, let me mention a debate 
that has been going on among DMRG practitioners on 
whether calling the DMRG algorithm a renormalization- 
group method is a misnomer or not. The answer depends 
on what one considers the essence of a RG method. If 
this essence is the systematic thinning out of degrees of 
freedom leading to effective Hamiltonians, DMRG is an 
RG method. However, it does not involve an ultraviolet 
or infrared energy cutoff in the degrees of freedom which 
is at the heart of traditional RG methods and hence often 
considered as part of the core concepts. 



A. Real-space renormalization of Hamiltonians 

The set of concepts grouped under the heading of 
"renormalization" has proven extremely powerful in pro- 
viding succinct descriptions of the collective behavior of 
systems of size L with a diverging number of degrees of 
freedom N^^^^. Starting from some microscopic Hamil- 
tonian, degrees of freedom are iteratively integrated out 
and accounted for by modifying the original Hamilto- 
nian. The new Hamiltonian will exhibit modified as well 
as new couplings, and renormalization group approxi- 
mations typically consist in physically motivated trun- 
cations of the set of couplings newly generated by the 
a priori exact elimination of degrees of freedom. One 
obtains a simplified ( "renormalized" ) effective Hamilto- 
nian that should catch the essential physics of the sys- 
tem under study. The success of this approach rests on 
scale separation: for continuous phase transitions, the di- 
verging correlation length sets a natural long-wavelength 
low-energy scale which dominates the physical proper- 
ties, and fluctuations on shorter length scales may be 
integrated out and summed up into quantitative modifi- 
cations of the long-wavelength behavior. In the Kondo 
problem, the width of the Kondo resonance sets an en- 
ergy scale such that the exponentially decaying contri- 
butions of energy levels far from the resonance can be 
integrated out. T his is the essence of the numerical renor- 
malization group llWilsonL [19751. 

Following up on the Kondo problem, it was hoped that 
thermodynamic-limit ground-state properties of other 
many-body problems such as the one-dimensional Hub- 
bard or Heisenberg models might be treated similarly, 
with lattice sites replacing energy levels. However, results 
of real-space renormalization schemes turned out to be 
poor. While a precise analysis of this observation is hard 
for many-body problems, ^hitc and Noack (1992 ) iden- 
tified the breakdown of the real-space renormalization 
group for the toy model of a single non-interacting parti- 
cle hopping on a discrete one-dimensional lattice ("par- 
ticle in a box" problem). For a box of size L, the Hilbert 
space spanned by {|j)} is L-dimensional; in state the 



compound block AA 



block A 



block A 



FIG. 1 Lowest-lying eigenstates of blocks A (dashed) and 
AA (solid) for the single particle in a box problem in the 
continuum limit. 



particle is on site i. The matrix elements of the Hamilto- 
nian are band-diagonal, {i\H\i) = 2; {i\H\i ± 1) = — 1 in 
some units. Consider now the following real-space renor- 
malization group procedure: 

1. Interactions on an initial sublattice ("block") A of 
length i are described by a block Hamiltonian Ha 
acting on an M-dimensional Hilbert space. 

2. Form a compound block A A of length 2£ and the 
Hamiltonian Haa, consisting of two block Hamil- 
tonians and interblock interactions. Haa has di- 
mension M^. 

3. Diagonalize Haa to find the M lowest-lying eigen- 
states. 

4. Project Haa onto the truncated space spanned by 
the M lowest-lying eigenstates, Haa 



rrtr 



5. Restart from step 2, with doubled block size: 2£ 
£, AA — > A, and Haa ~^ ^a, until the box size is 
reached. 



The key point is that the decimation procedure of the 
Hilbert space is to take the lowest-lying eigenstates of 
the compound block A A. This amounts to the assump- 
tion that the ground state of the entire box will essen- 
tially be composed of energetically low-lying states living 
on smaller blocks. The outlined real-space renormaliza- 
tion procedure gives very poor results. The breakdown 
can best be understood visually (Fig. ^ : assuming an al- 
ready rather large block size, where discretization can be 
neglected, the lowest-lying states of A all have nodes at 
the lattice ends, such that all product states of AA have 
nodes at the compound block center. The true ground 
state of AA has its maximum amplitude right there, such 
that it cannot be properly approximated by a restricted 
number of block states. 

Merely considering is olated blocks imposes wrong 
boundary conditions, and lWhite and Noackl l)l992l) could 
obtain excellent results by combining Hilbert spaces from 
low-lying states of block A obtained by assuming vari- 
ous combinations of fixed and free boundary conditions 
(i.e. enforcing a vanishing wave function or a vanishing 
wave function derivative respectively at the boundaries). 
They also realized that combining various boundary con- 
ditions for a single particle would translate to accounting 
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for fluctuations between blocks in the case of many in- 
teracting particles. This observation in mind, we now 
return to the original question of a many-body problem 
in the thermodynamic limit and formulate the following 
strategy: To analyze which states have to be retained for 
a finite-size block A, A has to be embedded in some en- 
vironment, mimicking the thermodynamic limit system 
A is ultimately embedded in. 



B. Density matrices and DIVIRG truncation 

Consider, instead of the exponentially fast growth 
procedure outlined above, the following linear growth 
prescription ijWhitel ^992): Assume that for a sys- 
tem (a block in DMRG language) of length i we 
have an A/'^-dimensional Hilbert space with states 
{|m|^)}. The Hamiltonian Hi is given by matrix elements 
{mf\Hi\rhf). Similarly we know the matrix representa- 
tions of local operators such as (mf |Cj |mf ). 

For linear growth, we now construct He+i in the prod- 
uct basis {Imfcr)} = {\m,f)\a^)}, where \a^) are the 
A^site local states of a new site added. 

The thermodynamic limit is now mimicked by embed- 
ding the system in an environment of the same size, as- 
sumed to have been constructed in analogy to the system. 
We thus arrive at a superblock of length 2^-1-2 (Fig. I^J, 
where the arrangement chosen is typical, but not manda- 
tory. 

As the final goal is the ground state in the thermo- 
dynamic limit, one studies the best approximation to it, 
the ground state of the superblock, obtained by numeri- 
cal diagonalization: 



IV') — ^ ] ^ ] ^ ^ ^ ^ 4'm° cr'^m 



(4) 



blocks block E 

\ 2 sites \ 

]o o[ 



superblock 



system 



environment 



T 



new block S 



new block E 



FIG. 2 System meets environment: Fundamental DMRG 
construction of a superblock from two blocks and two single 
sites. 



through a reduced density-matrix p. 



(5) 



where the states of the environment have been traced 
out, 



(6) 



p has N^' eigenvalues Wa and orthonormal eigenstates 
p\wa) = Wa\wa), with '^^Wa = 1 and Wa > 0. Wc 
assume the states are ordered such that wi > W2 > W3 > 
.... The intuition that the ground state of the system is 
best described by retaining those states with largest 
weight Wa in the density matrix can be formalized as 
follows. Consider some bounded operator A acting on 
the system, such as the energy per lattice bond; ||A|| = 
max0 |((/)|A|0)/(0|0)| = CA- The expectation value of A 
is found to be, using Eq. Q and Eq. ©, 



(A) - (V'|i|7/^)/(^|^) =Tr5/5A 



(7) 



Expressing Eq. ||7J) in the density-matrix eigenbasis, one 
finds 



where tp^s^s^E^E = {m^a^ia^m^ltp). {\m^a^)} = 
{\i)} and {\m^a^)} = {\j)} are the orthonormal prod- 
uct bases of system and environment (subscripts have 
been dropped) with dimensions — AI^'NgHc and 
= M^Nsitc respectively (for later generalizations, 
we allow ^ N^). Some truncation procedure from 
to < states must now be implemented. Let 
me present three lines of argumentation on the optimiza- 
tion of some quantum mechanical quantity, all leading 
to the same truncation prescription focused on density- 
matrix properties. This is to highlight different aspects 
of the DMRG algorithm and to give confidence in the 

prescription found. 

Optimization of expectation values ijWhitd . Il998fl : If 
the superblock is in a pure state as in Eq. Q), sta- 
tistical physics describes the physical state of the system 



(A) = ^ Wa{Wa\A\Wa 



(8) 



Then, if we project the system state space down to the 
AI^ dominant eigenvectors \wa) with the largest eigen- 
values. 



(^)approx = Wa{Wa\A\ 
a=l 

and the error for (A) is bounded by 



Wo 



(9) 



|(^)approx - (i)| < I E Wa\cA= ipCA- (10) 
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This estimate holds in particular for energies. Several re- 
marks are in order. Technically, I have neglected to trace 
the fate of the denominator in Eq. upon projection; 
the ensuing correction of (1 — Cp)^^ is of no relevance to 
the argument here as Cp —>■ 0. The estimate could be 
tightened for any specific operator, knowing {wa\A\wa) , 
and a more efficient truncation procedure be named. For 
arbitrary bounded operators acting on the system, the 
prescription to retain the state spanned by the dom- 
inant eigenstates is optimal. For local quantities, such as 
energy, magnetization or density, errors are of the order 
of the truncated weight 



(11) 



which emerges as the key estimate. Hence, a fast de- 
cay of density matrix eigenvalues Wa is essential for the 
performance of this truncation procedure. The trunca- 
tion error of Eq. (|10|) is the total error only if the system 
had been embedded in the final and exactly described 
environment. Considering the iterative system and envi- 
ronment growth and their approximate representation at 
each step, additional sources of an environmental error 
have to be considered. In practice, therefore, errors for 
observables calculated by DMRG are often much larger 
than the truncated weight, even after additional steps to 
eliminate environmental errors have been taken. Care- 
ful extrapolation of results in AI^ (rather ep) is therefore 
highly recommended, as will be disc ussed below. 

Optimization of the wave function jWhitilliillliil: 
Quantum mechanical objects are completely described 
by their wave function. It is thus a reasonable demand 
for a truncation procedure that the approximative wave 
function \ip) where the system space has been truncated 
to be spanned by only orthonormal states \a) = 



IV') = EE I") I-?')' 

Q=l j=l 

minimizes the distance in the quadratic norm 



(12) 



(13) 



This problem finds a very compact solution in a singu- 
lar value decomposition (SVD), which was the original 
approach of White; SVD will be considered in a slightly 
different setting in the next section. The following is an 
alternative, more pedestrian approach. Assuming real 
coefficients for simplicity, one has to minimize 



1 - 2 ^ tpijOajUai + ^ ( 



■a] 



(14) 



must be stationary for the global minimum of the dis- 
tance (I13f) . where we have introduced the density- matrix 
coefficients 



= E^y■^^'J■ 



(16) 



Equation (|15|l is stationary, according to the Rayleigh- 
Ritz principle, for \a) being the eigenvectors of the den- 
sity matrix. Expressing Eq. (|15|l in the density-matrix 
eigenbasis, the global minimum is given by choosing \a) 
to be the eigenvectors \wa) to the largest eigenvalues 
Wa of the density matrix, as they are all non-negative, 
and the minimal distance squared is, using Eq. Hill) . 



(17) 



a=l 



The truncation prescription now appears as a variational 

principle for the wave function. 

Optimization, of entanglement jCaite . J200l|, l2003t 
iLatorre et all 120041: lOsborne and Nielser , l2002|) : Con- 
sider the superblock state \tp) as in Eq. The essen- 
tial feature of a nonclassical state is its entanglement, 
the fact that it cannot be written as a simple prod- 
uct of one system and one environment state. Bipar- 
tite entanglement as relevant here can best be studied 
by repre senting in its form afte r a Schmidt decom- 
position ijNielsen and Chuanell200Ci|) : Assuming without 
loss of generality > N^, consider the (N^ x N^)- 
dimensional matrix A with Aij = tpij. Singular value 
decomposition guarantees A = UDV^, where U is 
{N^ X Af-^)-dimensional with orthonormal columns, D 
is a (N^ X Af-^)-dimensional diagonal matrix with non- 
negative entries Daa — \/Wa, and V'^ is a {N^ x A^-^)- 
dimensional unitary matrix; lip) can be written as 

JV-B jv-^ 

1^) = EEE^-v^^<S-N)b') (18) 

i—l a—1 j — 1 



TV* 



The orthonormality properties of U and ensure that 
kf> = J2iUia\i) and = Y^j^jalj) form orthonor- 

mal bases of system and environment respectively, in 
which the Schmidt decomposition 



1^) 



-/Vschmidt 

E 



/^hf)kf) 



(19) 



with respect to Oaj and Uai ■ To be stationary in Oaj , we 



must have ip. 



finding that 



1-E 



(15) 



holds. N^N^ coefficients are reduced to A^schmidt < 
non-zero coefficients ^Wq, wi > W2 > > .... 
Relaxing the assumption A^"^ > A^^, one has 



A^schmidt < min(A^^,Ar^). 



(20) 



7 



The suggestive labelling of states and coefficients in Eq. 
(|19|l is motivated by the observation that upon tracing 
out environment or system the reduced density matrices 
for system and environment are found to be 

Afschmidt -'Vschmidt 

a CK 

(21) 

Even if system and environment are different fSec. lII.D|l . 
both density matrices would have the same number of 
non-zero eigenvalues, bounded by the smaller of the di- 
mensions of system and environment, and an identical 
eigenvalue spectrum. Quantum information theory now 
provides a measure of entanglement through the von Neu- 
mann entropy, 

-^ Schm idt 

S'vN = -Tr/51n2 /5 = - ^ Wa\n2Wa- (22) 

a = l 

Hence, our truncation procedure to states preserves 
a maximum of system-environment entanglement if we 
retain the first AI^ states \wa) for a = 1, . . . , as 
— a;ln2a; grows monotonically for < x < 1/e, which 
is larger than typical discarded eigenvalues. Let me 
add that this optimization statement holds strictly only 
for the unnormalized truncated state. Truncation leads 
to a wave function norm ^1 — ep and enforces a new 
normalization which changes the retained Wq. and Sypf. 
While one may easily construct density-matrix spectra 
for which upon normalization the truncated state pro- 
duced by DMRG does no longer maximize entanglement, 
for typical DMRG density-matrix spectra the optimiza- 
tion statement still holds in practice. 

C. Infinite-system DIVIRG 

Collecting the results of the last two sections, the so- 
called in finite- system D MRG algorithm can now be for- 
mulated iWhitel Il992|) . Details of efficient implementa- 
tion will be discussed in subsequent sections; here, we 
assume that we are looking for the ground state: 

1. Consider a lattice of some small size i, forming the 
system block S. S lives on a Hilbert space of size 

with states {|M/)}; the Hamiltonian and 
the operators acting on the block are assumed to 
be known in this basis. At initialization, this may 
still be an exact basis of the block {N^^^^ < M^). 
Similarly, form an environment block E. 

2. Form a tentative new system block S' from S and 
one added site (Fig.|2J). S' lives on a Hilbert space of 
size = M^Nsitc, with a basis of product states 
{|M/(t)} = {|M/)|ct)}. In principle, the Hamil- 
tonian H£_^_i acting on S' can now be expressed in 
this basis (which will not be done explicitly for effi- 
ciency, see Sec. lH.I|) . A new environment E' is built 
from E in the same way. 



3. Build the superblock of length 2£ + 2 from S' and 
E'. The Hilbert space is of size N^'N^^, and the 
matrix elements of the Hamiltonian H21+2 could 
in principle be constructed explicitly, but this is 
avoided for efficiency reasons. 

4. Find by large sparse-matrix diagonalization of 
H21+2 the ground state {ip). This is the most time- 
consuming part of the algorithm fSec. IH.II1 . 

5. Form the reduced density- matrix p — Tr E\'ip){'fp\ as 
in Eq. © and determine its eigenbasis \wa) ordered 
by descending eigenvalues (weight) Wa ■ Form a new 
(reduced) basis for S' by taking the M^' eigenstates 
with the largest weights. In the product basis of 
S', their matrix elements are {mf a\m^^-^); taken as 
column vectors, they form a x M^' rectangular 
matrix T. Proceed likewise for the environment. 

6. Carry out the reduced basis transformation = 
T^H^+iT onto the new M'^-state basis and take 
Hl'i^i Hi+i for the system. Do the same for the 
environment and restart with step (2) with block 
size £+1 until some desired final length is reached. 
Operator representations also have to be updated 
(see Sec HTTTl . 

7. Calculate desired ground state properties (energies 
and correlators) from this step can also be car- 
ried out at each intermediate length. 

If the Hamiltonian is reflection-symmetric, one may con- 
sider system and environment to be identical. One is not 
restricted to choosing the ground state for lip); any state 
accessible by large sparse-matrix diagonalization of the 
superblock is allowed. Currently available algorithms for 
this diagonalization limit us however to the lowest-lying 
excitations (see Sec. III.ip . 



D. Infinite-system and finite-system DMRG 

For many problems, infinite-system DMRG does not 
yield satisfactory answers: The idea of simulating the fi- 
nal system size cannot be implemented well by a small 
environment block in the early DMRG steps. DMRG is 
usually canonical, working at fixed particle numbers for 
a given system size (Sec.^^J. Electronic systems where 
the particle number is growing during system growth to 
maintain particle density approximately constant are af- 
fected by a lack of "thermalization" of the particles in- 
jected during system growth; t-J models with a rela- 
tively small hole density or Hubbard models far from 
half-filling or with complicated filling factors are particu- 
larly affected. The strong physical effects of impurities or 
randomness in the Hamiltonian cannot be accounted for 
properly by infinite-system DMRG as the total Hamilto- 
nian is not yet known at intermediate steps. In systems 
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FIG. 3 Finite-system DMRG algorithm. Block growth and 
shrinkage. 



with strong magnetic fields or close to a first order tran- 
sition one may be trapped in a metastable state favored 
for small system sizes e.g. by edge effects. 

Finite-system DMRG manages to eliminate these con- 
cerns to a very large degree and to reduce the error (al- 
most) to the truncation error. The idea of the finite- 
system algorithm is to stop the infinite-system algorithm 
at some preselected superblock length L which is kept 
fixed. In subsequent DMRG steps (Fig. one applies 
the steps of infinite-system DMRG, but instead of simul- 
taneous growth of both blocks, growth of one block is 
accompanied by shrinkage of the other block. Reduced 
basis transformations are carried out only for the grow- 
ing block. Let the environment block grow at the expense 
of the system block; to describe it, system blocks of all 
sizes and operators acting on this block, expressed in the 
basis of that block, must have been stored previously 
(infinite-system stage or previous applications of finite- 
system DMRG). When the system block reaches some 
minimum size and becomes exact, growth direction is re- 
versed. The system block now grows at the expense of 
the environment. All basis states are chosen while system 
and environment arc embedded in the final system and 
in the knowledge of the full Hamiltonian. If the system 
is symmetric under refiection, blocks can be mirrored at 
equal size, otherwise the environment block is shrunk to 
some minimum and then regrown. A complete shrinkage 
and growth sequence for both blocks is called a sweep. 

One sweep takes about two (if reflection symmetry 
holds) or four times the CPU time of the starting infinite- 
system DMRG calculation. For better performance, 
there is a simple, but powerful "prediction algorithm" 
(see Sec. 111. Ill , which cuts down calculation times in finite- 
system DMRG by more than an order of magnitude. In 
fact, it will be seen that the speed-up is the larger the 
closer the current (ground) state is to the final, fully con- 
verged result. In practice, one therefore starts running 
the infinite-system DMRG with a rather small number of 
states Mq, increasing it while running through the sweeps 
to some final Mfinai S> Mq . The resulting slowing down of 
DMRG will be offset by the increasing performance of the 
prediction algorithm. While there is no guarantee that 
finite-system DMRG is not trapped in some metastable 
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state, it usually finds the best approximation to the 
ground state and convergence is gauged by comparing re- 
sults from sweep to sweep until they stablize. This may 
take from a few to several dozen sweeps, with electronic 
problems at incommensurate fillings and random poten- 
tial problems needing most. In rare cases it has been 
observed that seemingly converged finite-system results 
are again suddenly improved after some further sweeps 
without visible effect have been carried out, showing a 
metastable trapping. It is therefore advisable to carry 
out additional sweeps and to judge DMRG convergence 
by carrying out runs for various M. Where possible, 
choosing a clever sequence of finite-system DMRG steps 
may greatly improve the wave function. This has been 
successfully attempted in momentum space DMRG (Sec. 
IVII.A|) and quantum chemistry DMRG CSec. IVII.Bjl . 

To show the power of finite-system DMRG, consider 
the calculation of the plaquette currents induced on a t- 
J ladde r by imposing a source c urrent on one edge of the 
ladder l)Schollwock et all l2003|) . Fig. shows how the 
plaquette currents along the ladder evolve from sweep to 
sweep. While they are perfectly converged after 6 to 8 
sweeps, the final result, a modulated exponential decay, 
is far from what DMRG suggests after the first sweep, 
let alone in the infinite-system algorithm (not shown). 

Despite the general reliability of the finite-system al- 
gorithm there might be (relatively rare) situations where 
its results may be misleading. It may for example be 
doubted whether competing or coexisting types of long- 
range order are well described by DMRG, as we will 
see that it produces a very specific kind of wave func- 
tions, so-called matrix-product states fSec. lIII.X|l . These 
states show either long-range order or, more typically, 
short-ranged correlations. In the case of competing forms 
of long-range order, the infinite-system algorithm might 
preselect one of them incorrectly e.g. due to edge ef- 
fects, and the finite-system algorithm would then be quite 
likely to fail to "tunnel" to the other, correct kind of 
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long-range order due to the local nature of the improve- 
ments to the wave function. Perhaps such a problem is at 
the origin of th e current disagreernent between state-of- 
the-art DMRG fJcckelmann' 2002b. '2003b HZhangL 1200^ 
and QMC l,Sandvik et al... ,2004, 2003) on the extent of 
a bond-ordered wave phase between a spin-density-wave 
phase and a charge-density-wave phase in the half-filled 
extended Hubbard model; but no consensus has emerged 
yet. 



E. Symmetries and good quantum numbers 

A big advantage of DMRG is that the algorithm con- 
serves a variety of, but not all, symmetries and good 
quantum numbers of the Hamiltonians. They may be 
exploited to reduce storage and computation time and 
to thin out Hilbert space by decomposing it into a sum 
of sectors. DMRG is optimal for studying the lowest- 
lying states of such sectors, and refining any such de- 
composition immediately gives access to more low-lying 
states. Hence the use of symmetries is standard practise 
in DMRG implementations. Symmetries used in DMRG 
fall into three categories, continuous abelian, continuous 
nonabelian, and discrete. 

Continuous abelian symmetries. The most frequently 
implemented symmetries in DMRG are the C/(l) symme- 
tries leading to total magnetization S^^^ and total par- 
ticle number TVtot sls good (conserved) quantum num- 
bers. If present for some Hamiltonian, all operators can 
be expressed in matrix form as dense blocks of non-zero 
matrix elements with all other matrix elements zero. 
These blocks can be labelled by good quantum num- 
bers. In DMRG, reduced basis transformations preserve 
these block structures if one fixes total magnetization 
and/or total particle number for the superblock. As- 
sume that block and site states can at a given DMRG 
step be labeled by a good quantum number, say, particle 
number N. This is an essential prerequisite (cf. trans- 
lational invariance leading to momentum conservation; 
see below). As the total number of particles is fixed, we 
have iVtot = N^s + N^s + N„e + N,^e. Equation ® 
implies that only matrix elements {m^ a^\p\rh^ a^) with 
A'^^s -I- N^s = N^s + N^s can be non-zero. The density 
matrix thus has block structure, and its eigenvectors from 
which the next block's eigenbasis is formed can again be 
labeled by particle number N,-^s + N^s. Thus, for all 
operators only dense blocks of non-zero matrix elements 
have to be stored and considered. For superblock wave 
functions, a tensorial structure emerges as total particle 
number and/or magnetization dictate that only product 
states with compatible quantum numbers (i.e. adding up 
to the fixed total) may have non-zero coefficients. 

The performance gains from implementing just the 
simple additive quantum numbers magnetization and 
particle number are impressive. Both in memory and 
CPU time typically significantly more than an order of 
magnitude will be gained. 



Continuous nonabelian symmetries. Nonabelian sym- 
metries that have been considered are the qu antum group 
symmetry SUJ2) llsierra and Nishind.ll997l) SU( 2 ) spin 
symm e try fMcCulloch and Gulacsil l2000l I2OOII l2002t 
IWadaL I2OOQ: .Xiang et aL 2001), and the charge 5C/{2) 
pseudospin svmmetrv |McCulloch and GiilacsiL |2002), 
which holds for the bipartite Hubbard model without 
field (Yang and Zhang, 1990): its generators are given 

by 1+ = ^ = ^ud = 

Implementation of nonabelian symmetries is much 
more complicated than that of abelian symmetries , 
the most performant one llMcCulloch and Gulacsill2002l) 
building on Clcbsch-Gordan transformations and elimi- 
nation of quantum numbers via the Wigner-Eckart the- 
orem. It might be crucial for obtaining high-quality re- 
sults in applications with problematically large truncated 
weight such as in two dimensions, where the truncated 
weight is cut by several orders of magnitude compared 
to implementations using abelian symmetries only; the 
additional increase in performance is comparable to that 
due to use of [/(I) symmetries compared to using no sym- 
metries at all. 

Discrete symmetries. I shall formulate them for a 
fermionic Hamiltonian of spin-^ particles. 

Spin-flip symmetry: if the Hamiltonian H is invariant 
under a general spin flip t^^i, one may introduce the spin 
fiip operator P — Yi^ Pi, which is implemented locally on 
an electronic site as A|0) = |0), P,| t) = U), P^U) = | T), 
Pi\ U) = -\ U) (fermionic sign). As [H,P] = 0, there 
arc common eigenstates of H and P with P\ip) = ± | ''/')■ 

Particle-hole symmetry: if H is invariant under a 
particle-hole transformation, [H, J] = for the particle- 
hole operator J = Jli*^*' where the local operation is 
given by J, |0) = |n), J.|T) = (-1)^1 i), Mi) = ("1)1 T), 
Ji\ Ti) = |0)i introducing a distinction between odd and 
even site sublattices, and eigenvalues J\ip) = ±\ip). 

Reflection symmetry (parity): in the case of reflection 
symmetric Hamiltonians with open boundary conditions, 
parity is a good quantum number. The spatial reflec- 
tion symmetry operator C2 acts globally. Its action on a 
DMRG product state is given by 

C-alm^o-^CT^m^) = {-iy^\m^c7^a'^m^) (23) 

with a fermionic phase determined by 77 = {N^s -\- 
N„s){N^E -\- N^e). Again, eigenvalues are given by 
C2IV') — i I ''/')■ Parity is not a good quantum number 
accessible to finite-system DMRG, except for the DMRG 
step with identical system and environment. 

All three symmetries commute, and an arbitrary nor- 
malized wave function can be decomposed into eigen- 
states for any desired combination of eigenvalues ±1 by 
successively calculating 
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\^)±o\^) 



(24) 
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where d = P,J,C2. 

Parity may be easily implemented by starting the su- 
perblock diagonalization from a trial state fSec. lII.ljl that 
has been made (anti)symmetric under reflection by Eq. 
(I24|l . As [H,C2] = 0, the final eigenstate will have the 
same reflection symmetry. Spin-flip and particle-hole 
symmetries may be implemented in the same fashion, 
provided P and J are generated by DMRG; as they are 
products of local operators, this can be done along the 
lines of Sec. 111. (11 Another way of implementing these 
two local symmetries is to realize that the argument given 
for magnetization and particle number as good density- 
matrix quantum numbers carries over to the spin-flip and 
particle-hole eigenvalues such that they can also be im- 
plemented as block labels. 

Missing symmetries. Momentum is not a good quan- 
tum number in real-space DMRG, even if periodic bound- 
ary conditions (Sec. I11.H(I are used and translational in- 
variance holds. This is because the allowed discrete mo- 
menta change during the growth process and, more im- 
portantly, because momentum docs not exist as a good 
quantum number at the block level. Other DMRG vari- 
ants with momentum as a good quantum number will be 
considered in Sec. IVII.AI 

F. Energies: ground states and excitations 

As a method working in a subspace of the full Hilbert 
space, DMRG is variational in energy. It provides upper 
bounds for energies E{M) that improve monotonically 
with M, the number of basis states in the reduced Hilbert 
space. Two sources of errors have been identified, the en- 
vironmental error due to inadequate environment blocks, 
which can be amended using the finite-system DMRG 
algorithm, and the truncation error. Assuming that the 
environmental error (which is hard to quantify theoreti- 
cally) has been eliminated, i.e. finite-system DMRG has 
reached convergence after sufficient sweeping, the trun- 
cation error remains to be analyzed. Rerunning the cal- 
culation of a system of size L for various M, one observes 
for sufficiently large values of M that to a good approx- 
imation the error in energy per site scales linearly with 
the truncated weight, 

{E{M) ~ E,^,,^)/L cx ep, (25) 

with a non-universal proportionality factor typically of 
or der 1 to 10, somet i mes more (this observation goes back 
to lWhite and Husd l)l993(l : iLeeeza and FathI |)1996D give 
a careful analysis). As tp is often of order 10~^° or less, 
DMRG energies can thus be extrapolated using Eq. H25() 
quite reliably to the exact M — oo result, often almost at 
machine precision. The precision desired imposes the size 
of M, which for spin problems is typically in the lower 
hundreds, for electronic problems in the upper hundreds, 
for two-dimensional and momentum-space problems in 
the lower thousands. As an example for DMRG precision, 
consider the results obtained for the ground-state energy 



TABLE I Ground state energies per site Eq of the 5=1 
isotropic antiferromagnetic Heisenberg chain for M block 
state s kept and associated tr uncated weight ep(Af). Adapted 
from I White and Husd J1993D . 



M 


Eo{M) 


e,{M) 


36 


-1.40148379810 


5.61 X 10"* 


72 


-1.40148403632 


3.42 X 10"^° 


110 


-1.40148403887 


1.27 X 10"" 


180 


-1.401484038970 


1.4 X 10"^^ 



per site of the S = 1 antiferromagnetic Heisenberg chain 
in Table n 

Experiments relate to energy differences or relative or- 
dering of levels. This raises the question of calculating 
excitations in DMRG. Excitations are easiest to calculate 
if they are the ground state in some other symmetry sec- 
tor of the Hamiltonian and are thus algorithmically not 
different from true ground states. If, however, we are in- 
terested in some higher-lying state in a particular Hilbert 
space sector, DMRG restricts us to the lowest-lying such 
states because of the restrictions of large sparse-matrix 
diagonahzations (Sec. III.I|I . Excited states have to be 
"targeted" in the same way as the ground state. This 
means that they have to be calculated for the superblock 
at each iteration and to be represented optimally, i.e. re- 
duced basis states have to be chosen such that the error 
in the approximation is minimized. It can be shown quite 
easily that this amounts to considering the eigenstates of 
the reduced density-matrix 

Ps = Tr£^a,|V'.)(V'.|, (26) 

i 

where the sum runs over all targeted states (ground 
and a few excited states) and '^i — 1- There is no 
known optimal choice for the on, but it seems empirically 
to be most reasonable to weigh states roughly equally. To 
maintain a good overall description for all targeted states 
at a fixed M, typically less than 5 or so excited states are 
targeted. Best results are of course obtained by running 
DMRG for each energy level separately. 

G. Operators and correlations 

In general, we will also be interested in evaluating 
static n-point correlators O-"^ = O^^^ . . . with re- 
spect to some eigenstate of the Hamiltonian. The most 
relevant cases are n = 1 for density or local magnetiza- 
tion, and n = 2 for two-point density-density, spin-spin 
or creation-annihilation correlators, {niUj), {S^ S^) or 

(4s>. 

Let us first consider the case n — 1. The iterative 
growth strategy of DMRG imposes a natural three-step 
procedure of initializing, updating and evaluating corre- 
lators. 
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1. Initialization: Oi acts on site i. When site i is added 
to a block of length £ — 1, {a\Oi\a) is evaluated. With 
{Im^)} being the reduced basis of the new block incorpo- 
rating site i and {|m^_i)} that of the old block, one has 

{me\di\rhi) = ^ {■me\mi^ia)(a\d^\a){mi^ia\mi). 

(27) 

{niilnii-icr) is already known from the density-matrix 
eigenstates. 

2. Update: At each further DMRG step, an approxi- 
mate basis transformation for the block containing the 
site where Oi acts from {|mf)} to {jm^+i)} occurs. As 
Oi does not act on the new site, the operator transforms 
as 

{mi+i\d\mi+i) = 

{mi+i\mia) {me\d\mi) {rh£a\mi+i) . (28) 

This expression is evaluated efficiently by splitting it into 
two O(M^) matrix-matrix multiplications. 

3. Evaluation: After the last DMRG step 
{m^ (T^ m-^ \^p) is known and {Oi) reads, assum- 
ing Oi to act on some site in the system block, 

(v-io.i^) = 

(m^a^a^m^lV'). (29) 

In the case of 2-point correlators two cases have to 
be distinguished, whether the locations i and j of the 
contributing 1-point operators act on different blocks or 
on the same block at the last step. This expression is 
again evaluated efficiently by splitting it into two 0{M^) 
matrix-matrix multiplications. 

If they act on different blocks, one follows through 
the procedure for 1-point operators, yielding {m^\Oi\rh^) 
and (m^|Oj The evaluation is done by the follow- 
ing modification of the 1-point case, 

imo.i^p) = 

{■m^\d,\m^){m^a^a^m^\^). (30) 

If they act on the same block, it is wrong to obtain 
(to I Oi Oj I m) through 

(to|C),;Oj|to) = ^(to|6,|to')(to'|6j|to) (false), (31) 

where an approximate partition of unity has been in- 
serted. However, it is only one to a very good approxima- 
tion when it is used to project the targeted wave function, 
for which it was constructed, but not in general. 



Instead, such operators have to be built as a compound 
object at the moment when they live in a product Hilbert 
space, namely when one of the operators acts on a block 
(of length I — 1), the other on a single site, that is being 
attached to the block. Then we know {mi-i\Oi\rhi-i) 
and {a\0.j\a) and within the reduced bases of the block 
of length £ 

{mt\didj\frit) = Y, 

{me\rni-ia){m£-i\Oi\mi^i) x (32) 
{a\dj\a) {rhi_ia\rhi) 

is exact. Updating and final evaluation for "compound" 
operators proceed as for a one-point operator. 

One-point operators show similar convergence behav- 
ior in M as local energy, but at reduced precision. 

While there is no exact variational principle for two- 
point correlations, derived correlation lengths are mono- 
tonically increasing in M, but always underestimated. 
The underestimation can actually be quite severe and be 
of the order of several percent while the ground state en- 
ergy has already converged almost to machine precision. 

As DMRG always generates wave functions with ex- 
ponentially decaying correlations fSec. IIII..^ . pow er-law 
decay s of correlations are problematic. fXndcrsson et al\ 
(|1999) show that for free fermions the resulting corre- 
lation function mimics the correct power-law on short 
length scales increasing with M, but is purely exponen- 
tial on larger scales. However, the derived correlation 
length diverges roughly as M^-^, such that for M ^ oo 
criticality is recovered. 

H. Boundary Conditions 

From a physical point of view periodic boundary con- 
ditions are normally highly preferable to the open bound- 
ary conditions used so far for studying bulk properties, 
as surface effects are eliminated and finite-size extrapola- 
tion works for much smaller system sizes. In particular, 
open boundaries introduce charge or magnetization oscil- 
lations not always easily distingui shable from true cha rge 
density waves or dimerization fsee lWhite et all l)2002() for 
a thorough discussion on using bosonization to make the 
distinction). 

However, it has been observed early in the history of 
DMRG that ground state energies for a given M are much 
less precise in the case of periodic boundary conditions 
than for open boundary conditions with differences in 
the relative errors of up to several orders of magnitude. 
This is reflected in the spectrum of the reduced density- 
matrix, that decays much mor e slowly (see Sec . IIII.BII . 
However, it has been shown bv lVerstraete et aTl ()2004b|) 
that this is an artefact of the conventional DMRG setup 
and that, at some algorithmic cost, essentially the same 
precision for a given M can be achieved for periodic as 
for open boundary conditions fSec. IIII.'X|l . 
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FIG. 5 Typical system and environment growth for periodic 
boundary conditions in the infinite-system algorithm. 



To implement periodic boundary conditions in the 
infinite-system DMRG algorithm the block-site structure 
is typically changed as shown in Fig. [SJ other setups are 
however also feasible and used. For finite-system DMRG, 
the environment block grows at the expense of the sys- 
tem block, then the system block grows back, till the 
configuration of the end of the infinite-system algorithm 
is reached. This is repeated with changed roles (unless 
translational invariance allows to identify system and en- 
vironment at equal size) . A minor technical complication 
arises from the fact that blocks grow at both ends at var- 
ious steps of the algorithm. 

Beyond the usual advantages of periodic boundary con- 
ditions, combining results for periodic and antiperiodic 
boundary conditions allows the calculation of responses 
to boundary conditions such as spin stiffness, phase sen- 
sitivity iSchmitteckcrt and Eckern, 1996J or superfiuid 
density ijRapsch et a/.lll999f) . Periodic and antiperiodic 



boundary conditions c}^_|_]^ = ±c| are special cases of 

the general complex boundary condition c^_|_]^ = e"^c|. 
Implementing the latter is a tedious, but straightfor- 
ward generalization of real-valued DMRG; memory dou- 
bles, computation time quadruples. Numerical stabil- 
ity is assured because the density matrix remains her- 
mitian. This generalization has been used on a ring 
with interactions and impurities to determine the cur- 
rent /(</)) which is neither saw tooth-like nor sinusoidal 
ijMeden and SchoUwockL l2003a|) , to obtain the conduc- 
tance of interacting nanowires ijMeden and Schollwockl 
|2003b*). For open boundary conditions, complex- valued 
DMRG has been used to introduce infinitesimal current 
source terms for time-rev ersal symmetry breaking in elec- 
tronic ladder structures ijSchollwock et a/.l. 12005 



I. Large sparse matrix diagonalization 

Algorithms. Key to DMRG performance is 
the efficient diagonalization of the large sparse su- 
perblock Hamiltonian. All large sparse matrix di- 
agonalization algorithms iteratively calculate the de- 
sired eigenstate from some (random) starting state 



through successive costly matrix-vector multiplications. 
In DMRG, the two algorithms ty pically use d are 
the Lanczos method iCulhim and W illonghh^ 119851 
iGolub and van Loan. 1996) and the Jacobi-Davidson 
method (Sleiipen and van der Vorstl Il99l . The pleas- 
ant feature of these algorithms is that for a x di- 
mensional matrix it takes only a much smaller number 
A^ <C A of iterations such that iterative approximations 
to eigenvalues converge very rapidly to the maximum and 
minimum eigenvalues of H at machine precision. With 
slightly more effort other eigenvalues at the edge of the 
spectrum can also be computed. Typical values for the 
number of iterations (matrix-vector multiplications) in 
DMRG calculations are of the order of 100. 

Representation of the Hamiltonian. Naively, the su- 
perblock Hamiltonian is a M^A"j?j^g dimensional matrix. 
As matrix- vector multiplications scale as (dimension)^, 
DMRG would seem to be an algorithm of order 0{M*). 
In reality, it is only 0{M^), as typical tight-binding 
Hamiltonians act as sums over two operator terms: 
Assuming nearest-neighbor interactions, the superblock 
Hamiltonian decomposes as 



H — Hs + Hst + H,, + H,i 



Hp 



(33) 



Hs and He contain all interactions within the system 
and environment blocks respectively, are hence of dimen- 
sion M. Multiplying them to some state is of or- 
der M''Aj?;jg. Hs, and H,e contain interactions between 
blocks and the neighboring sites, hence are of dimension 
MNsitc- Consider a typical interaction S^S^_^-^, where £ 
is the last site of the block and ^ -I- 1 a single site. Then 

(34) 

and multiplying this term to {ip) is best carried out in a 
two step sequence: the expression 



(TOfjrn 



J2 {ma\S+S-\m'a'){m'a'Tn\ij), (35) 



that is of order 0{M^ N^^^^) for the determination of all 
state coefficients is decomposed as 

{m'aTn\iy) ^^{cr\S-\a'){m'a'Tn\iP), (36) 



of order ©(M^AtS^J and 



[m\S^\m') {m' (jTn\v) 



(37) 



of order 0{M^N'^:^^^), where an order of A'sitc is saved, im- 
portant for large Agite. The Hamiltonian is never explic- 
itly constructed. Such a decomposition is crucial when 
block-block interactions Hse appear for longer-ranged 
interactions. Considering again S'^ S~ ^ with S~ now act- 
ing on the environment block, factorization of the Hamil- 
tonian again allows to decompose the original term 

(77icrrn I (/)) — {mn\S^ \m' n') {m' am' \'ip) ^ (38) 
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that is of inconveniently high order 0{M'^N'^:^^^) for the 
determination of all state coefhcients, as 

{m'aTn\v) = ^(n|S'"|n')(r7i'cr™'|V'), (39) 

of order 0{APnI^^) and 

(m(Tm|0) — '^^{m\S~^\m'){m' aTn\v) ^ (40) 



of order 0{M^N'^^^^), such that an order of M is saved 
and the algorithm is generally of order 0{M^). This is 
essential for its performance. 

Eigenstate prediction. Providing a good guess for the 
final eigenstate as starting state of the iterative diago- 
nalization allows for arbitrary cutdowns in the number 
of iterations in the iterative diagonalization procedures. 

For both the infinite-system and the finite-system 
DMRG it is possible to provide starting states which of- 
ten have overlap w 1 with the final state, leading to a 
dramatic reduction of iterations often down to less than 
10, speeding up the algorithm by about an order of mag- 
nitude. 

In infinite-system DMRG, the physical system changes 
from step to step. It seems intuitive that for a very 
long system, the composition of the ground state from 
block and site states may only weakly depend on its 
length, such that the ground state coefficients remain 
almost the same under system growth. One might there- 
fore simply use the old ground state as prediction state. 
This fails; while the absolute values of coefficients hardly 
change from step to step in long systems, the block ba- 
sis states are fixed by the density-matrix diagonalization 
only up to the sign, such that the signs of the coeffi- 
cients are effectively random. Variou s ways of fixing these 

" )in and Louri200lt 



random sie 


ns have been proposed 


l^choUwocJ 


119981 


Sun e,t all 


|2002t) 



In the case of the finite-system DMRG, the physical 
system does not change from DMRG step to DMRG step, 
just th e struct ure of the effective Hilbert space changes. 
Iwhitd llmmh has given a prescription how to predict 
the ground state expressed in the block-site structure of 
the next DMRG step. If basis transformations were not 
incomplete, one could simply transform the ground state 
from one basis to the next to obtain a prediction state. 
The idea is to do this even though the transformation is 
incomplete. The state obtained turns out to be an often 
excellent approximation to the true ground state. 

Let us assume that we have a system with open bound- 
ary conditions, with a current system block of length £ 
and an environment block of length L — ^ — 2. The target 
state is known through {m£(7i-\-iae^2'niL-e-2\i') ■ Now as- 
sume that the system block is growing at the expense of 
the environment block, and we wish to predict the coeffi- 
cients (m£+itTf+2O'^+3TOL-£-3|'0)) where in both DMRG 
steps 1-0) is describing the same quantum mechanical 
state. As we also know {miai^i\me^i) from the current 
DMRG iteration and (mL_£_3cr£+3|mL_£_2) from some 



previous DMRG iteration, this allows to carry out two 
incomplete basis transformations: First, we transform 
system block and site to the new system block, giving 

{■mi+iai+2mL-e-2W = (41) 
{■me+i\meai+i){miae+iai+2mL-i-2\'ip); 

second, we expand the current environment block into a 
product state representation of a block with one site less 
and a site: 



{mi+iai+2cr£+3mL-£-3\^) = (42) 
{■mL-i-30-i+3\mL-i-2){mi+iai+2'mL-i-2\'4')- 



The calculation involves two 0{M^) matrix multipli- 
cations, and works independently of any assumptions for 
the underlying model. It relies on the fact that by con- 
struction our incomplete basis transformations are al- 
most exact when applied to the target state(s). 



J. Applications 

In this section, I want to give a very brief overview 
over applications of standard DMRG, which will not be 
covered in more detail in later sections. 

From the beginning, there has been a strong focus 
on one-dimensional Heisenberg models, in particular the 
5 = 1 case, where the Haldane g ap A = 0.41052 J was de - 
termined to 5 digit precision bv lWhite and Hus3 l|l993t ). 
Other autho rs have cons i dered t he equal time struc- 
ture fac tor ijSieling et all [2000; So rensen and Affleck . 
| l994albr) and focused on top ological order ijLou et al , 
l2003l I2OO2I: iQin et all \2Q0?i\ . There has been a par- 
ticular emphasis on the stud y of the existence of free 
S = ^ end spiri s in su c h chains llBatista et ajj.ll998lll99i: ; 



i Hallberg et^M. Il999l: riannod et all l200nt iPolizzi al 
mm IWhite and Husei [122^, where the open bound- 
ary conditions of DMRG are very useful; for the study 
of bulk properties authors typically attach real S — ^ 
end spins that bind to singlets with the effective ones, 
removing degeneracies due to boundary effects. Soon, 
studies carried over to the S = 2 Heisenberg chain, 
where the first reliable determination of the gap A = 
0.085(5) and the correlation length g » 50 was pro- 
vided by SchoUwock and Jolicoeur ('1995'), and confirmed 
and e n hanced in other works ( Aschauer and SchoUwock , 
1998t ICapone and Capraral I2OOII: iQin a/.l Il997t : 



SchoUwock et a;.Vl996aVSchollwock and JohcceurV199e 
Wada, 2000; Wang et al., 1999). The behavior of Hal- 
dane (integer sp in) chains in (stagger e d) ma g netic fields 
was studie d bv ISorense n and Affieckl ^ll993^■ I Lou et all 
(1999), Erc olessi et al\ l)2000(l . and lCapone and Caprarg 
(jSoOlJ). In general, DMRG has been very useful in the 
study of plateaux in magnetization processes of spin 
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chains 


llCitro et al'.. 


'2000: 'Hida. '200^: 'Kawaa;uchi et al. 


2002 




Lou et al.. 


.2000a.h.: .Silva- Valencia and Miranda 


2002 


:lTRndor, p.t oiJ. Il QPflHYa.mamnto F.t oil 1200(11). 



The isotropic half-integer spin Heisenberg chains are 
critical. The logarithmic corrections to the power- 
law of spin-spin correl ations in the S = ^ chain 

were first cons idered by 'Hallberg et al! (1995') , later by 



and important generalization of the Heisenberg model 
DMRG h as allowed full clarifi cation of the rich phas e 
diagram llAmmo n and I madal 12 001: It oi et all l200i 
Pati and Singh'.'2000 HPati et ad . ll998HYamashita et al 
2b00a. 1998, 2000b.). 



i 



Hikiha,ra, a,nd Fnrusa.k i (1998). ,Tsai and Marston (.2000}. 
Shiroishi et al\ l|200ll)~ and iBoos Rtal\ (|2002^ . For the 

S' = I case, the central charge c = 1 and again the loga- 
rithmic corrections to the sp in-spi n correlations were de- 



termined bv lHallberg et al 



chains were considered by 
case of transverse fields by 
Bilinear-biquadratic S 



199' 



I 111999 



. Qua s iperi odic S = 



Hieida et all l)2001 ) 



20001). and the 



DMRG in its finite-system version has also been very 
successful in the study of systems with impurities or ran- 
domness, where the true ground state can be found with 
reasonable precision . Impu ri ties have been studied by 
[Wang and Mallwitz! ('1996'), 'Schmitteckert and Eckern] 
11996), Mikcsk a et al (1997) Martins et al (1997), 
Schmitteckert e t'al\ Il99^, U^ aukamp et al\ (11998), 
Ng et aL (.2000(1 . iLou et a/.l ((1998,) . .Zhang et al\ (|l997^) . 
Zhan^t g/J (|l998b|) . whereas other authors have focused 



1 spi n 

extensive ly stud ie d llBursill et all . Il99 
ISchoUwoc k et 

iO, Il996bl) as well as the effect of 
Dzyaloshinskii-Moriya interactions l)Zhao et all 1200,^) . 

Important experimentally relevant generalizations of 
the Heisenberg model are obtained by adding frustrat- 
ing interactions or dimerization, the latter modelling 
either static lattice distortions or phonons in the adia- 



chains have been on random interactions or fields (iHida, 1996, 1997a 
ISatd. Il998t 



kolezhuk et al! 


("1996", 


'l997); 


Kolezhuk and SchoUwocH 


(2002 


); IMaeshima and Okunishi 


(2000); 


Pati et oL 


(199e 
(199e 


ll997al):IUhrig et all (ll999alhD: IWhite and Affieckl 

) have extensively studied the ground state phase 



diagrams of such systems. DMRG has been instrumental 
in the discovery and description of gap ped and gapless 
chiral ph ases in frust r ated sp i n chains lHikiliara|j2001 
12002; Hi kihara et Hi] . l2000a[ l2001albl : iKaburagi eHil 
Il999,) . Critical exponents for a supersymmetric spin 
chain were obtained by Scnthil et al 1 1999i) . 

As a first step towards two dimensions and due 
to the many experimental realizations, spin lad- 
ders have been among the first DMRG applica- 
tions goin g beyond si m ple H eisenberg chains, start- 
ing with IWhite et al\ (|1994() . In the meantime, 
DMRG has e merged as standa r d tool for th e study of 
spin ladders (jPd tl^^ [2001: Hikihara and Furusakl 
' 200lHKaw aguchi et ad 120 03^ .Legeza and Solvom. .199^ 
Trumper an d Gazza, 200 ll IWand. l200d: IWang et al\ . 
2002t IWhitel . il996at IZhu et all l200ll) . A focus 
of recent interest has been the effec t of cyclic ex- 
change interactions on spi n lad ders llHikihara et all 
12003; Honda and Horiguchi", "2001; 'LauchU et all \2003t. 
iNunner e t _aL .2,01)2: Sakai an d Hasegawa, J.999) • 

spin systerns. ferrimagnets have bee n 



Among other 



(T997b*c), 



'Tonegawa et al} lll99i 



studied by Pati et al. 
iLangari et al. (2000), and Hikihara et al. (2000b^. One- 
dimensional to y models of the Kagom e lattice have been 
investigated by 'Pati and Singh' 1*1999^). * White and Singhl 
(|2000|) . and AValdtraann et al ( 2OOO0, whereas su- 
persymmetric spin cha ins have been considered by 
iMarston and Tsail l|l999fl . 

Spin-orbit chains with spin and pseudospin degrees 
of freedom are the large-JJ limit of the two-band Hub- 
bard model at quarter-filling and are hence an interesting 



1999at iHikihara et aZ.L Il999t lJuozapavicius et al.- .1999 
Urba and RosengrenL l2003l) . There has also been 
work on edges and impu rities in Luttinger liq- 
uids llBediirftig et all \199A iQin et all Il996i Il997at 
ISchonhammer a/.l.l2000^. 

The study of electronic models is somewhat more com- 
plicated because of the larger number of degrees of free- 
dom and because of the fermionic sign. However, DMRG 
is free of the negative sign problem of quantum Monte 
Carlo and hence the method of choice for one-dimensiona l 
Hubbard c hains llAebischer et al 



electronic models. 




Afigia et al!, '200(J!; 'DauJ, '200 (^_ lDaylj,n d Noack. 
1998, 2000 ; Daul a nd Scalapin 

ILepetit et all I2OO0I : iLepetit and Pastoi 



, IP 

iMaurel and Lepetid . l2000t iNishimoto et al 
ISakamoto and KubnL 1199(1 IZhausi Il997ri^d 



Hubbard lad ders ^o nca et al, 2000 t iDaul et 
2000b; Hamacher et ah. 2002: Marston et all 12002 : 



Noack et all Il997l Il994 Il995albl: IScalapino et al 



200^: 'SchoDwock et al\ 1200.4 IVoita et all Il999l 12001 



Weihong et al., .2001) have been studied in the en- 
tire range of interactions as the precision of DMRG 
is found to depend only very weakly on the inter- 
action strength. The three-band Hub b ard m odel 
has been considered by iJeckelmann et al\ l|l998|) and 
INishimoto et al\ (|2002b|) . Similarl y, authors have stud- 
ied b o th the t-J model on c hains (IChen and Moukouri 



19961 iDoublet and Lepctit. 


"1999' 'Maurel et al!, |200ll 


Mutou et al.. 


.1998a: 


White and Scalapino 


.1997bj) and 


on ladders jHavwood et all Il995l: iRommer et a/.l. I200C 


Siller et all l200ll l2002l IWhite and Sralapinrl n997a 



1998c^ 



The persistent current response to m agnetic fluxes 
throug h rin gs has be en studied by 'Bv rnes et al\ 
(|2002() and iMeden and SchoUwock (2003aj, also in 
view of possible quasi-exact conductance calcula tions 
(|Meden and SchollwockL l2003bl: iMoIina et all^OO^ . 

Even before the advent of the Hubbard model, 
a very similar model, the Pariser- Parr-Pople model 
(iPariser and Pari! Il953t iPopleL [l95l, accommodating 
both dimerization effects and longer-range Coulomb 
interaction, had been invented in quantum chemistry 
to study conjugated polymer structures. DMRG 
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auth 

mi 

199c 


ors llAnusoova et all 
; IBarford et all 
: iBoman and Bursill 


Jl997| 
|2» 

19 


IBarford and Bursill 
[ |Benda,zzoli et all 

98; Kuwabara et al;, 


1998 


: iLeoetit and Pastor. 


.1997. 
200(1 




Raghu et al.,. 


■2002t 


Ramasesha et all Il997l 






I19974 


19981 Il997hld: lYaron et a.l 


J. 11998 


) , moving to 


more 



and more reali stic models of_ molecules, such as 
Dolvdiacatvl ene ijRace et aZI . I2 OOI, 20 0^), polyfpara- 
phcnylene) ('Anusoo^^e^^j 

[l997i .Barford et all Il998t 

Bursill and Barford, 2002), or pol y enes \'QsxiaYde^il. , 
"oOlt IBursill and BarfordL Il999t IZhang and George , 
200 1|) . Polymer phonon s have been considere d also in 
the non-adiabatic case l|Barford et all l2002a|) . There 
is closely related work on the Pei erls-Hubbard model 
llAnusoova et ad ll999tlJeckelmannlil998; Qtsuka. ,199& 
iPang and LiangLll995(l 

Sin ce the early days of DMRG history l|Yu and WhiteL 
Il993() interest 
tice, 
erant 
both 
199 
200 



has also focused 
generic one-dimensional 
electrons and localized 
for the one-channel 



JCarruzo_an^^ , 
iMcCulloch et al 




on the Kondo lat- 
structures of itin- 
magnetic moments, 
Tara and RoscngrejL 



1997albL Il996bl Il997ct Isikkema eit oi.L Il99 
1998t IWa tanabe etall Il999fl and two-channel case 



Garcia all _ I2OOC 
IShibata et al^ 



1996a. 



Wane . 



llMorenq ' et a/.,. ( 20011). An d erson m odels have been stud 
ied b y IGuerrero and Yul lll99c 1 . IGuerrero and Noackl 



l|1996 f). and G uerrero and Noackl l|200H) . 

The boso nic version of the Hu b bard m odel has been 
studie d by iKiihner and MonienI ||1998|) . iKtihner et al\ 
l|2000l) . and'KoUath et al} ('2004'), with emphasis on dis- 
order effects by Rapsch et al. (1999). 

Going beyond one dimension, two-dimensional 
electron gases in a Landau level have been mapped 
to one-dimensional m odels suitable for D MRG 
jBcrgholt z and Karlhed'A l2003t IShibata and YoshiokaL 
I2OOL .pool: DMRG was a pplied to molecular iron 
rings l|Normand et all l200l|) and has elucidated the 
lowest rota tional band of the g iant K eplerate molecule 
MoraPeso ()Exler and SchnackL |2003|) . More generic 
higher-dimensional applications will be discussed later. 

Revisiting its origins, DMRG can also be used to 
provide high-ac curacy solutions in one-part icle quan- 
tum mechanics l)Martm-Delgado et all Il999|) : as there 
is no entanglement in the one-particle wave function, 
the reduced basis transformation is formed from the 
lowest-lying states of the superblock projected onto the 
system (after reorthonormalization) . It has been ap- 
plied to an asymptotically free model in two dimensions 
l|Martm-Delgado an d Sierra. 1999) and mo dified for up 
to three dimensions (jMartfn-Delgado et all I2OOLI . 



III. DMRG THEORY 

DMRG practitioners usually adopt a quite pragmatic 
approach when applying DMRG to study some physical 



system. They consider the convergence of DMRG results 
under tuning the standard DMRG control parameters, 
system size L, size of the reduced block Hilbert space 
M , and the number of finite-system sweeps, and judge 
DMRG results to be reliable or not. Beyond empiricism, 
in recent years a coherent theoretical picture of the con- 
vergence properties and the algorithmic nature of DMRG 
has emerged, and it is fair to say that we have by now 
good foundations of a DMRG theory: DMRG generi- 
cally produces a particular kind of ansatz states, known 
in statistical physics as matrix-product states; if they 
well approximate the true state of the system, DMRG 
will perform well. In fact, it turns out to be rewarding 
to reformulate DMRG in terms of variational optimiza- 
tion within classes of matrix-product states fSec. IIII.'X|l . 
In practice, DMRG performance is best studied by con- 
sidering the decay of the eigenvalue spectrum of the re- 
duced density-matrix, which is fast for one-dimensional 
gapped quantum systems, but generically slow for criti- 
cal systems in one dimension and all systems in higher 
dimensions (Sec. 1111. B|) . This renders DMRG applica- 
tions in such situations delicate at best. A very coherent 
understanding of these properties is now emerging in the 
framework of bipartite entanglement measures in quan- 
tum information theory. 



A. Matrix-product states 

Like conventional RG methods, DMRG builds on 
Hilbert space decimation. There is however no Hamil- 
tonian flow to some flxed point, and no emergence of 
relevant and irrelevant operators. Instead, there is a flow 
to some fixed point in the space of the reduced density 
matrices. As has been pointed out by various authors 
llDukelskv et a/.Ul998tfMartfn-Delgado and SierraLll996l: 
Ostlund and Rommer, 1995; Romme r and OstlundL 
1997; Takasaki et al . J999) . this implies that DMRG 
enerates po s ition- d ependent matr i x-prod uct states 
Fannes et all . 119891: iKliimoer et'oB. Il993|) as block 
states. However, there are subtle, but crucial difi'er- 



ences between DMRG s tates and matrix- p roduct states 
l|Takasaki et all Il999t IVerstraete et all l2004bj) that 
have important consequences regarding the variational 
nature of DMRG. 

Matrix-product states are simple generalizations of 
product states of local states, which we take to be on 
a chain, 



(43) 



obtained by introducing linear operators ^^[(7^] depend- 
ing on the local state. These operators map from some 
A/-dimensional auxiliary state space spanned by an or- 
thonormal basis to another M-dimensional auxil- 

iary state space spanned by {|a)}: 



(44) 



a/3 
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One may visualize the auxiliary state spaces to be located 
on the bonds (i,i + 1) and (i — 1,?). The operators are 
thus represented by M x Af matrices {Ai[ai])ai3; M will 
be seen later to be the number of block states in DMRG. 
We further demand for reasons explained below that 



(45) 



A position-dependent unnormalized matrix-product 
state for a one-dimensional system of size L is then given 

by 

IV') = E ((-^iin^'^NI^fl)) k), (46) 

{or} \ i=l I 

where and \4>b) are left and right boundary states 
in the auxiliary state spaces located in the above visual- 
ization to the left of the first and to the right of the last 
site. They are used to obtain scalar coefficients. Position- 
independent matrix-product states are obtained by mak- 
ing Eq. H44|l position- independent, Ai\ui\ At\a^. For 
simplicity, we shall consider only those in the following. 

For periodic boundary conditions, boundary states are 
replaced by tracing the matrix-product: 



{cr} 



(47) 



The best-known matrix-product state is the valence- 
bond-solid ground state of the bilinear-biquadratic 
S = \ Af5eck-Kenn edy-Lieb-Tasaki Hamiltonian 
(lAfheck el a/J . 119871 Il98j^ . where M = 2. 

Correlations in matrix-produ ct states: Conside r two lo - 
cal bosonic (for simplicity; see lAndersson et oZI l)l999() ') 
operators Oj and Oj+z, acting on sites j and j + l, applied 
to the periodic boundary condition matrix-product state 
of Eq. g7I). The correlator C{1) = (V'lOjOj+i 
is then found, with TrXTrF = Ty{X Y) and (ABC) (g) 
{XYZ) = {A(g)X){B(g) Y){C ® Z), to be given by 



Trf 



(48) 



where we hav e used the following mapping 
l|Andersson et al\ . Il999l: |R ommer and Ostlundl Il997j) 

from an operator O acting on the local state space to 
a M^-dimensional operator O acting on products of 
auxihary states \a(3) = \a) (g) \f3): 



O = 



a'\d\(j)A*[a'](E)A[ 



(49) 



Note that A* stands for A complex-conjugated only as 
opposed to A^ . Evaluating Eq. H48I) in the eigenbasis of 
the mapped identity, I, we find that in the thermody- 
namic limit L — > cxD 



C{1) 



exp(-Z/6) 



(50) 



with = — l/ln|Ai|. The Ai are the eigenvalues of I, 
and the Ci depend on O. This expression holds because 
due to Eq. (|^ |A.j| < 1 and Ai = 1 for the eigen- 
state (q;/3|Ai) = dap- Equation is thus seen to en- 
sure normalizability of matrix product states in the ther- 
modynamic limit. Generally, all correlations in matrix- 
product states are either long-ranged or purely exponen- 
tially decaying. They are thus not suited to describe 
critical behavior in the thermodynamic limit. Even for 
gapped one-dimensional quantum systems their utility 
may seem limited as the correlators C{1) of these systems 
are generically of a two-dimensional classical Ornstein- 
Zernike form. 



C{1) 



Vi 



(51) 



whereas exponential decay as in Eq. 150|) is typical of one- 
dimensional classical Ornstein-Zcrnike forms. However, 
matrix-product states such as the AKLT ground state 
arise as quantum disorder points in general Hamiltonian 
spaces as the quantum remnants of classical phase tran- 
sitions in two-dimensional classical systems; these disor- 
der points are characterized by dimensional reduction of 
their correlations and typically characterize the qualita- 
tive properties of subsets of the Hamiltonian space, which 
turns them into most useful toy models, as e xemplified by 
the AKLT state l|Schollwock et al\ . \l996^i . Away from 
the disorder points, choosing increasingly large M as di- 
mension of the ansatz matrices allows to model the true 
correlation form as a superposition of exponentials for 
increasingly large I; even for power-law correlations, this 
modelling works for not too long distances. 

DMRG and matrix-product states: To show that a 
DMRG calculation retaining M block states produces 
M X AI matrix-product states, lOsthmd and R.ommed 
l)l995j) considered the reduced basis transformation to 
obtain the block of size i, 

(52) 

such that 

|to£_i) (g) Icr^.). (53) 

The reduced basis transformation matrices Ai[ai] auto- 
matically obey Eq. which here ensures that {|m^)} 
is an orthonormal basis provided {|m£_i)} is one, too. 
We may now use Eq. (|53|l for a backward recursion to 
express \me-i) via |m£_2) and so forth. There is a (con- 
ceptually irrelevant) complication as the number of block 
states for very short blocks is less than M. For simplicity, 
I assume that N^^.^ = AI, and stop the recursion at the 
shortest block of size N that has AI states, such that 
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where we have boundary-site states \mfj) = \ai . . . crff), 
hence 

ai,...,at 

(55) 

A comparison to Eq. 146|) shows that DMRG gener- 
ates position-dependent M x M matrix-product states 
as block states for a reduced Hilbert space of M states; 
the auxihary state space to a local state space is given 
by the Hilbert space of the block to which the local site 
is the latest attachment. Combining Eqs. 10} and (|55|l . 
the superblock ground state of the full chain is the vari- 
ational optimum in a space spanned by products of two 
local states and two matrix-product states, 

m^m^ {<J } 

{Al/2-i[^L/2-1\ ■ ■ ■ ^N+lW N+lVlmS ^(ai...af,) X 
(^L/2+2[ctl/2+2] • ■ • ^L-Jv[^L-Jv])m,-^,(^i+i_jv---'^i) 
|(Ti... 0-l), (56) 

which I have written for the case of the two single sites at 
the chain center; an analogous form holds for all stages 
of the finite-system algorithm, too. 

For gapped quantum systems we may assume that for 
blocks of length l':^ the reduced basis transformation 
becomes site-independent, such that the ansatz matrix A 
generated by DMRG is essentially position-independent 
in the bulk. At criticality, the finite-dimensional matrix- 
product state generated introduces some effective corre- 
lation length (growing with M). In fact, this has been 
numerically verified for free fermions at criticality, where 
the ansatz matrices for bulk sites converged exponentially 
fast to a position-independent ansatz matrix, but where 
this convergence slowed down with M l|Andersson et al\ . 
[l999). 

The effect of the finite-system algorithm can be seen 
from Eq. H5()l) to be a sequence of local optimization steps 
of the wave function that have two effects: on the one 
hand, the variational coefficients V^mSo-So-Bm-E are opti- 
mized, on the other hand, a new improved ansatz matrix 
is obtained for the growing block, using the improved 
variational coefficients for a new reduced basis transfor- 
mation. 

In practical applications one observes that even for 
translationally invariant systems with periodic bound- 
ary conditions and repeated applications of finite-system 
sweeps the position dependency of the matrix-product 
state does not go away completely as it str ictly should, in- 
dicati ng roo m for further im p rovem ent. iDukelskv et al\ 
l|l998l ) and 'Tak asaki etoR (|^9^ have pointed out 
and numerically demonstrated that finite-system DMRG 
(and TMRG; see Sec. IVIIl)) results can be improved and 
better matrix product states for translationally invari- 
ant Hamiltonians be produced by switching, after con- 
vergence is reached, from the S»»E scheme for the finite- 
system algorithm to a S»E scheme and to carry out 



some more sweeps. The rationale is that the variational 
ansatz of Eq. H56|l generates (after the Schmidt decom- 
position and before truncation) ansatz matrices of di- 
mension MN site at the two local sites due to Eq. I|2U|1 . 
whereas they are of dimension M at all other sites; this 
introduces a notable position dependence, deteriorating 
the overall wave function. In the new scheme, the new 
ansatz matrix, without truncation, has also dimension M 
(the dimension of the environment in Eq. (|20|l ) , such that 
the local state is not favored. The variational state now 
approaches its global optimum without further trunca- 
tion, just improvements of the ansatz matrices. 

The observation that DMRG produces variational 
states formed from products of local ansatz matrices has 
inspired the construction of variational ansatz states for 
the ground state of Hamiltonians (recurre nt variational 
approa ch; see Marti'n-Delgado and Sierra in lPeschel et al\ 
(p99a)) and the dominant eigenstate of transfer matri- 
ces (ten sor product variationa l appro ac h: .Gendiar et al 
mOS) : iGendiar and Nishind ll2002tl : iMaeshima et all 

il2nnih : lNishino et ai\ il2nniL l2nnfft). 

Closer in spirit to the original DMRG concept is the 
product wave function renormalization group or P WFRG 
l|Hieida et flIl . ll997tlNishino and Okunishl Il995l) . which 
has been applied successfully to the magnetization pro- 
cess of spin chains in external field where infinite- 
system DMRG is highly prone t o metastable trap- 
ping feieida et a l, _1997, ,2001; Ok unishi et all Il999at 
ISato and Akuts^j. ir99fil) and to the rest r icted solid-on- 
solid ni odel ijAkutsu and AkutsuL Il998t lAkutsu et all 
bOOlalhT ). While in DMRG the focus is to determine for 
each iteration (superblock size) the wave function nu- 
merically as precisely as possible in order to derive the 
reduced basis transformation, PWFRG directly operates 
on the matrices A itself: at each iteration, one starts 
with an approximation to the wave function and local 
ansatz matrices related to it by a Schmidt decomposi- 
tion. The wave function is then somewhat improved by 
carrying out a few Lanczos steps at moderate effort and 
Schmidt-decomposed again. The transformation matri- 
ces thus obtained are used to transform (improve) the 
local ansatz matrices, which in turn are combined with 
the decomposition weights to define a wave function for 
the next larger superblock. The final result is reached 
when both the ansatz and transformation matrices be- 
come identical as they should at the DMRG fixed point 
for reduced basis transformations. 

Variational optimization in matrix-product states. If 
we compare Eq. to Eq. lf5H|l , we see that the DMRG 
state is different from a true matrix-product state to de- 
scribe \ip): The A matrices link auxiliary state spaces 
of a bond to the right of a site to those on the left for 
sites in the left block, but vice versa in the right block. 
This may be mended by a transposition. This done, one 
may write the prefactors of |(Ti . . . (7^) as a true product 
of matrices by rewriting i'rnSaL/2C'L/2+im.'^ as a A/ x M 
matrix (5'[cri/2Ci/2-i-i])mS.mE. The remaining anomaly 
is, as pointed out above, that formal "translational in- 
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FIG. 6 Pictorial representation of S»E (above) and S»»E 
(below) DMR G schemes for open bou ndary conditions as 
introduced bv IVerstraete et al\ lj2004lJ l. (Paired) dots rep- 
resent auxiliary state spaces linked to a local state space. 
Straight lines symbolize maximal entanglement, ellipses and 
rectangles maps to local state spaces as detailed in the text. 
Note the special role of the boundary sites. Adapted from 
IVerstraete et al\ 



variance" of this state is broken by the indexing of 
by two sites, suggesting the modification of the S»»E 
scheme for the finite-system algorithm to a S«E scheme 
discussed above. However, as Vcrstrae te et 
have demonstrated, it is conceptually and algorithmically 
worthwhile to rephrase DMRG consistently in terms of 
matrix-product states right from the beginning, thereby 
also abandoning the block concept. 

To this end, they introduce (with the exception of the 
first and last site; see Fig. two auxiliary state spaces 
of dimension M, ai to the left and hi to the right of site 

such that on bond H. one has auxiliary state spaces bi 
and a^+i. They now consider maps from ag ^ bg ^ TCi 
from the product of two auxiliary state spaces to the 
local state space which can be written using matrices 
{A4a]U: A, = E., Ea,ft(^4^])a/3kf)(«^/3d- The 
|a) and \(3) are states of the auxiliary state spaces. 
On the first and last site, the corresponding maps map 
only from one auxiliary state space. These maps can 
now be used to generate a matrix p roduct state. To 
this purpose, IVerstraete et al\ l|2004bj) apply the string 
of maps Ai ® A2 ® ■ ■ ■ ® Al to the product of maxi- 
mally entangled states \4>i)\4>2) ■ ■ ■ where \cf)i) — 
E/3 ■ The maximal entanglement, in the 
language of matrix product states, ensures that the pref- 
actors of \ai . . .ai) are given by products of A[a'] ma- 
trices, hence this construction is a matrix-product state. 
Comparing this state to the representation of \il>) in Eq. 
(|56|1 . one finds that the maps A are identical to the (possi- 
bly transposed) basis transformation matrices A[a\ with 
the exception of the position of the single sites: in the 
S»»E setup (bottom half of Fig. EJ), there are no aux- 
iliary state spaces between the two single sites and one 
map corresponds to ^'[crfcrf+i]. In the S»E setup (top 
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FIG. 7 Periodic boundary condit ion setup used in the algo- 
rithm of IVerstraete et al\ i200 4br). Labelling as in previous 
figure; adapted from IVerstraete et al\ (120040) . 



half of Fig. ISJ, this anomaly disappears and correctly 
normalized A[a£] can be formed from ^'[cr^]. 

The (finite-system) DMRG algorithm for the S»E 
setup can now be reformulated in this picture as follows: 
sweeping forward and backward through the chain, one 
keeps for site £ all A at other sites fixed and seeks the 
"^i that minimizes the total energy. From this, one deter- 
mines Ai and moves to the next site, seeking V'^+i, and so 
on, until all matrices have converged. In the final matrix- 
product state one evaluates correlators as in Eq. H48|) . It 
is important to note that in this setup there is no trun- 
cation error, as explained above in the language of the 
Schmidt decomposition. Shifting the "active" site there- 
fore does not change the energy, and the next minimiza- 
tion can only decrease the energy (or keep it constant). 
This setup is therefore truly variational in the space of 
the states generated by the maps A and reaches a min- 
imum of energy within that space (there is of course no 
guarantee to reach the global minimum) . By comparison, 
the setup S»»E leads to a reduced basis transformation 
and always excludes two different auxiliary state spaces 
from the minimization procedure. It is hence not strictly 
variational. 

In this setup the generalization to periodic boundary 
conditions is now easy. Additional auxiliary state spaces 
ai and 5i are introduced and maximally entangled as for 
all other bonds. There is now complete formal transla- 
tional invariance of the ansatz (see Fig.]?)). On this setup, 
one optimizes maps (matrices) A one by one, going for- 
w ard and backwar d . 

IVerstraete et ll2n04bD have shown that for a given 
M, they obtain roughly the same precision for peri- 
odic boundary conditions as for open boundary condi- 
tions. This compares extremely favorably with standard 
DMRG for periodic boundary conditions where (in the 
worst case) up to states are needed for the same pre- 
cision. 



B. Properties of DMRG density matrices 

In order to gain a theoretical understanding of DMRG 
performance, we now take a look at the properties of 
the reduced density matrices and their truncation. Ob- 
viously, the ordered eigenvalue spectrum Wa of the re- 
duced density-matrix p should decay as quickly as possi- 
ble to minimize the truncated weight Cp = 1 — Eq=i 
for optimal DMRG performance. This intuitively clear 



19 



statement can be quantified: there are four major classes 
of density-matrix spectra, in descending order of DMRG 
performance. 

(i) Density-matrix spectra for M x M matrix-product 
states as exact eigenstates of quantum systems, with a 
finite fixed number of non- vanishing eigenvalues, leading 
to optimal DMRG performance. 

(ii) Density-matrix spectra for non-matrix-product 
states of one-dimensional quantum systems with expo- 
nentially decaying correlations, with leading exponential 
decay of Wa', spectra remain essentially unchanged for 
system sizes in excess of the correlation length. 

(iii) Density-matrix spectra for states of one- 
dimensional quantum systems at criticality, with a de- 
cay of Wa that slows down with increasing system size, 
leading to DMRG failure to obtain thermodynamic limit 
behavior. 

(iv) Density-matrix spectra for states of two- 
dimensional quantum systems both at and away from 
criticality, where the number of eigenvalues to be retained 
to keep a fixed truncation error grows exponentially with 
system size, restricting DMRG to very small system sizes. 

All scenarios translate to classical systems of one addi- 
tional dimension, due to the standard quantum-classical 
mapping from d- to {d + l)-dimensional systems. 

DMRG applied to matrix-product states. A state 
of the matrix-product form of Eq. (|46|l with dimension 
M, can be written as 

M M 
Q— 1 a—1 

where we have arbitrarily cut the chain into (left) system 
and (right) environment with 

E('^5in^NHk^>, (58) 

{o-s} ies 

and similarly |V'q); IV'f): and \^pa) a-''^ the cor- 

responding normalized states. An appropriate treat- 
ment of boundary sites is tacitly implied [cf. the discus- 
sion before Eq. H55I) ]. Then the density matrix ps — 

'}2^=i'*^a\i^a) {''l'a\ ^'^'^ ^as & finite spectrum of M non- 
vanishing eigenvalues. The truncated weight will thus be 
zero if we choose M > M for DMRG, as DMRG gener- 
ates these states [see Eq. 

In such cases, DMRG may be expected to become an 
exact method up to small numerical inaccuracies. This 
has been obs e rved re currently; Kaulke and Peschel (in 
iPeschel et ah ()l999a(l ) provide an excellent example, the 
non-hermitean g-symmetric Heisenberg model with an 
additional boundary term. This Hamiltonian is known 
to have matrix-product ground states of varying com- 
plexity (i.e. matri x sizes M) for partic ular choices of pa- 
rameters (see also lAlcaraz et a/l ()1994|) '). Monitoring the 
eigenvalue spectrum of the DMRG density matrix for a 
sufficiently long system, they found that indeed it col- 
lapses at these particular values: all eigenvalues but the 
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FIG. 8 Largest density-matrix eigenvalues vs. a tuning pa- 
rameter a in a. q = 1/4-symmetric Heisenberg chain. Largest 
eigenvalue « 1 invisible on logarithmic scale; eigenvalue col- 
lapse indicates pure matrix-product s tates of small finite d i- 
mension. From Kaulke and Peschel in lPeschel et al 
Reprinted with permission. 

M largest vanish at these points (Fig. |SJ). Similarly, a 
class of two-dimensional quantum Hamiltonians with ex- 
act matrix-produc t ground sta tes has been studied using 
DMRG e± a]\ ill 999ft . 

DMRG applied to generic gapped one- dimensional sys- 
tems. It is quite easy to observe numerically that for 
gapped one-dimensional quantum systems the eigenvalue 
spectrum of the density matrices decays essentially expo- 
nentially such that the truncated weight can be reduced 
exponentially fast by increasing M, which is the hallmark 
of DMRG success. Moreover, the eigenvalue spectrum 
c onverges to some therm odynamic-limit form. 

iPeschel et a,l\ (|l999lJ l have confirmed these numeri- 
cal observations by studying exact density-matrix spec- 
tra that may be calculated for the one-dimensional Ising 
model in a transverse field and the XXZ Heisenberg chain 
in its antifcrromagnetic gapped regim e, using a corner 
transfer matrix method llBaxteil[l98l . 

In those cases, the eigenvalues of p are given, up to a 
global normalization, as 

/ ^ \ 
w oc exp — E] Cj'^j (59) 
V j=o J 

with "fermionic" occupation numbers Uj = 0, 1 and an 
essentially linear "energy spectrum" e^: these are typi- 
cally some integer multiple of a fundamental scale e. The 
density-matrix eigenvalue spectrum shows clear exponen- 
tial decay only for large eigenvalues due to the increas- 
ing degree of degeneracy of the eigenvalue spectrum, as 
the number of possible partitions of en into ejUj grows. 
DMRG density matrices perfectly reproduce this behav- 
ior (Fig. 13 . Combining such exactly known corner trans- 
fer m atrix spectra with result s from the theory of parti- 
tions, lOkunishi et~ai\ l)l999bj) have derived the asymp- 
totic form 

Wa ^ exp (—const, x In^ a) (60) 
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FIG. 9 Density-matrix eigenvalues Wn vs. eigenvalue number 
n for a gapped XXY Heisenberg chain of L = 98 and three 
values of anisotropy A > Ac = 1. Degeneracie s are as pre- 
dicted analytically. From iPeschel et al\ (ll999b^ . Reprinted 
with permission. 

for the ath eigenvalue (see also lChan et 

These observations imply that for a desired truncated 
weight the number of states to be kept remains finite 
even in the thermodynamic limit and that the truncated 
weight decays exponentially fast in M , with some price 
to be paid due to the increasing degree of degeneracy 
occuring for large M. 

DMRG applied to one- dimensional systems at criti- 
cality. Numerically, one observes at criticality that the 
eigenvalue spectrum decays dramatically slower and that 
for increasin g system size this phen omenon tends to ag- 
gravate (jChung and Pesclieil l200l[l . with errors in e.g. 
ground-state energies increasing by sev eral orders of mag- 
nitude compared to gapped systems ijLegeza and FathL 
|l99S). 

Hence, the double question of the decay of the eigen- 
value spectrum for the density matrix of a given sys- 
tem size and the s i ze de pendency of this result arises. 
IChung and Peschell l|200l|) have shown for generic Hamil- 
tonians quadratic in fermionic operators that the density 
matrix spectrum is once again of the form of Eq. (|59|l . 
but the "energies" Cj are now no longer given by a simple 
relationship linear in j. Instead, they show much slower, 
curved growth, that slows down with system size. Trans- 
lating this into actual eigenvalues of the density matrix, 
they show less-than exponential decay slowing down with 
system size. This implies that for one-dimensional quan- 
tum systems at criticality numerical convergence for a 
fixed system size will no longer be exponentially fast in 
M. Maintaining a desired truncated weight in the ther- 
modynamic limit implies a diverging number M of states 
to be kept. 

DMRG in two-dimensional quantum systems. Due 
to the large interest in two-dimensional quantum sys- 



tems, we now turn to the question of DMRG conver- 
gence in g apped and critica l syste ms. In the early days 
of DMRG, iLiang and Pand l|l994r) had observed numer- 
ically that to maintain a given precision, an exponen- 
tially growing number of states M ~ a^, a > 1, had 
to be kept for system sizes L x L. However, reliable in- 
formation from numerics is very difhcult to obtain here, 
due to the very small syst em sizes in actual calculations. 
IChung and Peschell <)2000(l have studied a (gapped) sys- 
tem of interacting harmonic oscillators, where the density 
matrix can be written as the bosonic equivalent of Eq. 
(|59|l . with eigenvalues, again up to a normalization, 

w oc exp ^- ^ <^jb]b^ . (61) 

Considering strip systems of size L x N with N < L, 
numerical evaluations have shown that 

ej ~ const. X j/N, (62) 

hence the eigenvalue decay slows down exponentially 
with inverse system size. This 1/iV behavior can be 
understood by considering the spectrum of N chains 
without interchain interaction, which is given by the 
spectrum of the single chain with iV-fold degener- 
acy introduced. Interaction will lift this degeneracy, 
but not fundamentally change the slowdown this im- 
pose s on the decay of the density-ma trix spectrum [see 
also Idu Croo de Jongh and van Leeuwc n (1998) for the 
generic argument]. Taking bosonic combinatorics into 
account, one finds 

Wa ~ exp (— (const. /iV) In^ a) , (63) 

which is a consistent extension of Eq. (|50|l . For a sys- 
tem of size 10 X 10, typical truncation errors of 10^^ for 
M = 100, 10-^ for M = 500 and 10"^ for M = 1000 
have been reported (cf. Fig. I10|l . reflecting the very slow 
convergence of DMRG in this case. 

For a critical system, the non-interacting fermion 
model may be used once again as a model system 
fChun g and Peschel, 2001.) . For M = 2000 states, 
for this simple model, the resulting truncation error is 
5 X 10~^ for systems of size 12 x 12, 5 x 10~^ for 16 x 16 
and IQ-^ for size 20 x 20. Here, DMRG is clearly at its 
limits. 

DMRG precision for periodic boundary conditions. 
While for periodic boundary conditions the overall prop- 
erties of DMRG density matrices are the same as those 
of their open boundary condition counterparts, their 
spectra have been observed numerically to decay much 
more slowly. Away from criticality, this is due to some 
(usually only approximate) additional factor two degen- 
eracy of the eigenvalues. This can be explained by 
studying the amplitudes of density-matrix eigenstates. 
[Chung and Peschel ( 200(|) have demonstrated in their 
non-critical harmonic oscillator model that the density- 
matrix eigenstates associated to high eigenvalue weight 
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FIG. 10 Density-matrix eigenvalues ui„ for rectangular two- 
dimensional gapped systems of interacting harmonic oscilla- 
tors. The le ft-most curve corresponds to the one-dimensional 
case. From lChung and Peschell i2000l) . Reprinted with per- 
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FIG. 11 Amplitudes of highest-weight single-particle eigen- 
states of the left-block density-matrix for a one-dimensional 
open chain of L = 32 inte racting harmonic oscillators. From 
IChung and Peschell i2000D . Reprinted with permission. 



are strongly located close to the boundary between sys- 
tem and environment (Fig. Ill() . Hence, in a periodic 
system, where there are two boundary points, there are 
for the high-weight eigenvalues two essentially identical 
sets of eigenstates localized at the left and right boundary 
respectively, leading to the approximate double degener- 
acy of high-weight eigenvalues. At criticality, no such 
simple argument holds, but DMRG is similarly affected 
by a slower decay of spectra. 

In a recent study, IVerstraete et all l)2004b(l have shown 
that this strong deterioration of DMRG is essentially due 
to its particular setup for simulating periodic boundary 
conditions, and provided a new formulation of the algo- 
rithm which produces results of the same quality as for 
open boundary conditions for the same number of states 
kept, at the cost of losing matrix sparseness (see Sec. 

MM- 




FIG. 12 Diverging von Neumann entropy Sl vs. block length 
L for a critical isotropic XX chain in external fields H < 
Hcrit- Increasing fi eld strengths suppress entropy. From 
iLatorre et al\ (|22q3)- Reprinted with permission. 



C. DMRG and quantum information theory 

As under stood in th e early phases of DMRG 

development ijWhitel Il992t IWhite and NoackL [l99Z) , the 
reason for the success of the method is that no system 
is considered in isolation, but embedded in a larger en- 
tity. In fact, as discussed in Sec. III. B1 DMRG truncation 
can be understood in the language of quantum informa- 
tion theory as preserving the maximum entanglement be- 
tween system and environment as measured by the von 
Neumann entropy of entanglement, 



S = -Trp ln2 /5 = - ^ Wq ln2 Wa 



(64) 



in an M-dim ensional block state space. 

iLatorre ~e t al. ( 2004) have calculated the entropy of en- 
tanglement Sl for systems of length L embedded in infi- 
nite (an)isotropic XY chains and for systems embedded 
in finite, periodic Heisenberg chains of length iV, both 
with and without external field. Both models show criti- 
cal and noncritical behavior depending on anisotropy and 
field strength. In the limit N oo, Sl < i, which is 
obtained by observing that entropy is maximal if all 2^ 
states are equally occupied with amplitude 2^^. They 
find that Sl ^ oo for L ^ cx) at criticality, but satu- 
rates as Sl ^ S*j^ for L « ^ in the non-critical regime. 
At criticality, however, Sl grows far less than permitted 
Sl < i, but obeys logarithmic behavior 



Sl ^ k log2 L + const. 



(65) 



The constant k is found to be given by fc = 1/6 for the 
anisotropic XY model at the critical field He = 1, and 
by fc = 1/3 for the isotropic XY model a,t He = 1 as well 
as the isotropic Heisenberg model for H < He — 1 (for 
an isotropic XX model in field H < Hcrit, see Fig. I12|l . 
Away from criticality, the saturation value S^ decreases 
with decreasing correlation length ^ fFig. I13|l . 

One-dimensional quantum spin chains at criticality are 
described by conformally invariant ( 1 + 1 )-dime nsional 
field t heories. In fa c t, Eg . can be finked ijGaitel 
I2OO3I: ILatorre et~al[ \2004) to the geometric entropy 
l)Callan and Wilczekl . ll994() of such field theories, 



5r = ^iog2L, 



(66) 
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FIG. 13 Saturating density-matrix entropy Sl vs. block 
length L for an Ising chain in a transverse field H < Hait- 
Saturation entropy grows with H Hcth', divergence is re - 
covered at criticality (top curve). From lLatorre et al\ l|2004^ . 
Reprinted with permission. 

where c (c) are the central charges, if one observes that 
for the anisotropic XY model c = c — 1/2 and for the 
Heisenberg model and isotropic XY model c = c = 1. 

Geometric entropy arguments for {d + 1) -dimensional 
field theories use a bipartition of d-dimensional space by 
a (c? — l)-dimensional hypersurface, which is shared by 
system S and environment E. By the Schmidt decompo- 
sition, S and E share the same reduced density-matrix 
spectrum, hence entanglement entropy, which is now ar- 
gued to reside essentially on the shared hypersurface (cf. 
the locus of highest weight density-matrix eigenstates in 
Fig. ini see also|Gaitc (2001)). Taking the thermody- 
namic (infrared) limit, entropy scales as the hypersurface 
area, 



where A is some ultraviolet cutoff which in condensed 
matter physics we may fix at some lattice spacing. In- 
troducing a gap (mass), an essentially infrared property, 
into this field theory does not modify this behavior gener- 
ated on ultraviolet scales on the hypersurface. In d = 1, 
a more careful argument shows that there is a (sublead- 
ing) logarithmic correction as given above at criticality, 
saturation otherwise. 

Sl is the number of qubits corresponding to the en- 
tanglement information. To code this information in 
DMRG, one needs a system Hilbert space of size M > 
2Sl. jj-^ fg^Qi^ numerical results indicate that M and 2^'^ 
are - to a good approximation - proportional. Taking the 
linear dimensions of total space and embedded system to 
N,L oo, quantum information theory hence provides 
us with a unified picture of DMRG performance, which 
is in perfect agreement with the observations obtained by 
studying density matrix spectra: 

(i) In one-dimensional quantum systems away from 
criticality, DMRG yields very precise results for the ther- 
modynamic limit for some finite number of states kept, 
M r^2^L, which grows with the correlation length. 

(ii) In one-dimensional quantum systems at criticality, 
the number of states that has to be kept, will diverge as 

M{L) - L'', (68) 



with k from Eq. . This explains the failure of DMRG 
for critical one-dimensional systems as L — > oo. As fc is 
small, this statement has to be qualified; DMRG still 
works for rather large finite systems. 

(iii) In higher-dimensional quantum systems, the num- 
ber of states to be kept will diverge as 

M{L)^2^''~\ (69) 

rendering the understanding of thermodynamic limit be- 
havior by conventional DMRG quite impossible; infor- 
mation is beyond retrieval just as in a black hole - whose 
entropy scales with its surface, like the entanglement 
entropy would in a three-dimensional DMRG applica- 
tion! In any case, even for higher-dimensional systems, 
DMRG may be a very useful method as long as system 
size is kept resolutely finite, such as in nuclear physics 
or quantum chemistry applications. Recent proposals 
l|Verstraete and Cirad 1200 4') also show that it is possible 
to formulate generalized DMRG ansatz states in such a 
way that entropy shows correct size dependency in two- 
di mensional sys t ems (se e Sec. IVI|) . 

iLegeza et al\ ()2003a(l have carried the analysis of 
DMRG state selection using entanglement entropy even 
further, arguing that the acceptable truncated weight 
- which is not identical to, but closely related to 
the entropy of entanglement, which emerges as the 
key quantity ~ should be kept fixed during DMRG 
runs, determining how many states M have to be re- 
tained at each iteration. This dynamical block state 
select ion has already been applied in various contexts 
(Legeza e ^_aZjj20 03bt iLeeez a and Solvom, 2003). More 
recently, ^egezg^nd SolvomI (||200JD have tightened the 
relationship between quantum information theory and 
DMRG state selection by proposing a further refinement 
of state selection. M is now chosen variably to keep loss 
of quantum information below some acceptable thresh- 
old. They argue that this loss is given as 

X{S) = S{p) - PtypS{ptyp) - (1 - Ptyp)S{patyp)- (70) 

Here, p = ptypPtyp + (1 - Ptyp)Patyp- For a given M, 
Ptyp is formed from the M dominant eigenstates of p, 
Patyp from the remaining ones, with 1 — ptyp being the 
truncated weight and Pt,m .n.tvTi scaled to have trace one. 
iLeeeza and SolvomI l)2004|) report a very clear linear re- 
lationship between DMRG errors and x{£)- 

IV. ZERO-TEMPERATURE DYNAMICS 

As we have seen, DMRG is an excellent method to cal- 
culate ground states and selected excited eigenstates at 
almost machine precision. On the other hand, the target- 
ing of specific states suggests that DMRG is not suitable 
to calculate dynamical properties of strongly correlated 
systems even at T = as the time evolution of general ex- 
cited states will explore large parts of the Hilbert space. 
Closer inspection has revealed, however, that the rele- 
vant parts of Hilbert space can be properly addressed by 
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DMRG. For some operator A, one may define a (time- 
dependent) Green's function at T = in the Heisenberg 
picture by 



iGAit' -t)^{o\AHt')Ait)\o) 



(71) 



with t' > t for a time-independent Hamiltonian H. Going 
to frequency-range, the Green's function reads 



1 



Ea+LU + ir]-H 



A\0), (72) 



where rj is some positive number to be taken to zero at 
the end. We may also use the spectral or Lehmann rep- 
resentation of correlations in the eigenbasis of H, 

CaH = ^ |(n|i|0)|2(5(c^ + Eo- En). (73) 



This frequency-dependent correlation function is related 
to Ga {oj + iiy) as 

Ca(oj) = lim --ImGAii^ + iv)- (74) 

?;^0+ TT 

In the following, I shall also use Ga{(.^) and Ca{uj + irj) 
where the limit 77 — )■ 0^ will be assumed to have been 
taken in the former and omitted in the latter case. The 
role of 77 in DMRG calculations is threefold: First, it 
ensures causality in Eq. (|72|l . Second, it introduces a 
finite lifetime t cx I/77 to excitations, such that on finite 
systems they can be prevented from travelling to the open 
boundaries where their reflection would induce spurious 
effects. Third, 77 provides a Lorentzian broadening of 
Ca(c^), 



CAiuj + i7]) = ^ J dJGA^oj') 



V 



{U - LU')^ + jf 



(75) 



which serves either to broaden the numerically obtained 
discrete spectrum of finite systems into some "thermo- 
dynamic limit" behavior or to broaden analytical results 
for Ca for comparison to numerical spectra where 77 > 0. 

Most DMRG approaches to dynamical correlations 
center on the evaluation of Eq. H72|l . The first, which I 
ref er to as Lanczos vector dynamics, has been pioneered 
bv iHallberd l)l995j) . and calculates highly time-efficient, 
but comparatively rough approximations to dynami- 
cal quantities adopting the Balseiro-Gagliano method 
to DMRG. The second appr oach, which is referred to 
as correction vec t or me thod (|Kiihner and Whit3 . Il999t 
iRamasesha et al\ . llQQTj) . is yet another older scheme 
adopted to DMRG, much more precise, but numerically 
by far more expensive. A third approach, called DDMRG 
dynam ical DMRG), has been proposed bv Jcckclm amil 
2002aD : while on the surface it is only a minor reformu- 
lation of the correction vector method, it will be seen to 
be much more precise. 

Very recently, dynamical correlations as in Eq. H71|) 
have also been calculated directly in real time using 



time-dependent DMRG wi t h ada p tive Hilbert spaces 
llDalev et all 120041: IVidal l2004 IWhite and FeiguinL 
I2OO4D as will be discussed in Sec. IIX.GI These may then 
be Fourier transformed to a frequency representation. As 
time-dependent DMRG is best suited for short times (i.e. 
high frequencies) and the methods discussed in this Sec- 
tion are typically best for low frequencies, this alternative 
approach may be very attractive to cover wider frequency 
ranges. 

All approaches share the need for precision control and 
the elimination of finite-system and/or boundary effects. 
Beyond DMRG specific checks of convergence, precision 
may be checked by comparisons to independently ob- 
tained equal-time correlations using the following sum 
rules: 



(76) 
(77) 
(78) 



dujGA{Lo) = (0|iti|0) 
dujuuGAioj) = (0|i^[i/,i]|0) 
duuj^CAH = mA\H][H,A]\0), 



where the first equation holds for 77 > and the latter 
two only as 7; ^ 0. As DMRG is much more precise 
for equal-time than dynamical correlations, comparisons 
are made to DMRG results which for the purpose may 
be considered exact. Finite-size effects due to bound- 
ary conditions can be treated in various ways. They can 
be excluded completely by the use of periodic bound- 
ary co nditions at the price of much lower DMRG pre- 
cision llHallbersllT99l . In cases where open boundary 
conditions are preferred, two situations should be distin- 
guished. If A acts locally, such as in the calculation of an 
optical conductivity, one may exploit that finite ri expo- 
nentially suppresses excitations t,Jeckelmann. 20Q2a) . As 
they travel at some speed c through the system, a thermo- 
dynamic limit L — > 00, 77 — > with 77 = c/L may be taken 
consistently. For the calculation of dynamical structure 
functions such as obtained in elastic neutron scattering, 
A is a spatially delocalized Fourier transform, and an- 
other approach must be taken. The open boundaries in- 
troduce both genuine edge effects and a hard cut to the 
wave functions of excited states in real space, leading to 
a large spread in momentum space. To limit bandwidth 
in momentum space, filtering is necessary. The filtering 
function should be narrow in momentum space and broad 
in rea l space, while simultaneous ly strictly excluding edge 
sites. iKiihner and White! lll999D have introduced the so- 
called Parzen filter Fp , 



Fp(x) = 



1-6|2:|2-|-6|2:|3 0<|a;|<l/2 



2(1 -N)3 



1/2 < |x| < 1 



(79) 



where x — i/ Lp £ [—1,1], the relative site position in the 
filter for a total filter width 2Lp. In momentum space 
-Fp has a wave vector uncertainty Ag = 2-\/3/ip, which 
scales as if one scales Lp with system size L. For 



24 



finite-size extrapolations it has to be ensured that filter 
size does not introduce a new size dependency. This can 
be ensured by introducing a Parzen filter prefactor given 
by ^1407r/151ip. 



A. Continued fraction dynamics 

The technique of continued fraction dyn a mics has first 
been exploited by iGagliano and Balseird (^^3) in the 
framework of exact ground state diagonalization. Ob- 
viously, the calculation of Green's functions as in Eq. 
()72|l involves the inversion of H (or more precisely, 
Eq + Lu + it] — H), a typically very large sparse hermi- 
tian matrix. This inversion is carried out in two, at least 
formally, exact steps. First, an iterative basis transfor- 
mation taking _ff to a tridiagonal form is carried out. 
Second, this tridiagonal matrix is then inverted, allowing 
the evaluation of Eq. (|72|1 . 

Let us call the diagonal elements of H in the tridi- 
agonal form a„ and the subdiagonal elements 5^. The 
coefficients a„, 6^ are obtained as the Schmidt-Gram co- 
effcients in the generation of a Krylov subspace of unnor- 
malized states starting from some arbitrary state, which 
we take to be the excited state ^|0): 



with 



|/„+i>=i7|/„)-a„|/„)-fe'l/n-i>, 



l/o> = i|0), 

^ (/n|g|/n) 
{fn\fn) 

^2 ^ {fn-l\H\fn) _ 

" (/ri-l|/n-l) 



{fn\fn 



(80) 

(81) 
(82) 

(83) 



The global orthogonality of the states |/„) (at least in 
formal mathematics) and the tridiagonality of the new 
representation (i.e. {fi\H\fj) — for \i — j\ > 1) follow 
by induction. It can then be shown quite easily by an 
expansion of determinants that the inversion of _Eo + ^ + 
17] ~H leads to a continued fraction such that the Green's 
function Ga reads 



Ga{z) 



(0|iti|o) 



(84) 



where z = Eo+uj + ir]. This expression can now be evalu- 
ated numerically, giving access to dynamical correlations. 

In practice, several limitations occur. The iterative 
generation of the coefficients a„, &^ is equivalent to a 
Lanczos diagonalization of H with starting vector A\0). 
Typically, the convergence of the lowest eigenvalue of the 
transformed tridiagonal Hamiltonian to the ground state 
eigenvalue of H will have happened after n ^ O(IO^) 
iteration steps for standard model Hamiltonians. Lanc- 
zos convergence is however accompanied by numerical 



loss of global orthogonality which computationally is en- 
sured only locally, invalidating the inversion procedure. 
The generation of coeflnc i ents h as to be stopped before 
that. iKiihner and Whit'3 \199^ 

have proposed to mon- 
itor, for normalized vectors, (/o|/n) > e as termination 
criterion. The precision of this approach therefore de- 
pends on whether the continued fraction has sufficiently 
converged at termination. With A\0) as starting vector, 
convergence will be fast if A|0) is a long-lived excitation 
(close to an eigenstate) such as would be the case if the 
excitation is part of an excitation band; this will typically 
not be the case if it is part of an excitation continuuum. 

It should also be mentioned that beyond the above 
complications also arising for exact diagonalization with 
an exact H, additional approximations are introduced in 
DMRG as the Hamiltonian itself is of course not exact. 

Instead of evaluating the continued fraction of Eq. 
(I84II . one may also exploit that upon normalization of 
the Lanczos vectors |/„) and accompanying rescaling of 
the a„ and 6^, the Hamiltonian is iteratively transformed 
into a tridiagonal form in a new approximate orthonor- 
mal basis. Transforming the basis {|/n)} by a diago- 
nalization of the tridiagonal Hamiltonian matrix to the 
approximate energy eigenbasis of H, {\n)} with eigenen- 
ergies En, the Green's function can be written within this 
approximation as 



GAiuJ + irj) = 



(85) 



1 



En 



-UJ + b] - En 



|n)(n|i|0). 



where the sum runs over all approximate eigenstates. 
The dynamical correlation function is then given by 



Ca(w -I- it]) 



E 



\{n\m\ 



{Eo + uj - En)^ + 



(86) 



where the matrix elements in the numerator are simply 
the l/o) expansion coefficients of the approximate eigen- 
states \n). 

For a given effective Hilbert space dimension M, opti- 
mal precision within these constraints can be obtained by 
targeting not only the ground state, but also a selected 
number of states of the Krylov sequence. While a first 
approach is to take arbitrarily the first n states generated 
(say, 5 to 10) at equal weight, the approximate eigenbasis 
formulation gives direct access to the relative importance 
of the vectors. The importance of a Lanczos vector 
is given by, writing \A) = A\0), 



A„ = ^|(A|m)n(m|/„)|^ 



(87) 



which assesses the contribution the vector makes to an 
approximative eigenstate |to), weighted by that eigen- 
state's contribution to the Green's function. The density 
matrix may then be constructed from the ground state 
and a number of Lanczos vectors weighted according to 
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FIG. 14 Single magnon line of the S = 1 Heisenberg 
AFM from exact diagonalization, quantum Monte Carlo and 
DMR G for various system siz es and boundary conditions. 
From lKjilmer and Whitj il999tl . Reprinted with permission. 



FIG. 15 Spectral weight S+{q = tt,uj) of the 5 = 1/2 Heisen- 
berg AFM from Lanczos vector and correction vector DMRG. 
Nl indicates the number of target states; M = 256. Note that 
spectral weight times u is shown. From iKiihner and Whit3 
119930 • Reprinted with permission. 



Eq. (|87|) . The weight distribution indicates the apph- 
cability of the Lanczos vector approach: for the spin-1 
Heisenberg chain at q — tt, where there is a one-magnon 
band at a; = A = 0.41J, three target states carry 98.9 
% of total weight, for the spin- 1/2 Heisenberg chain at 
q = TT with a two-spinon continuum for < lu < ttJ, 
the first three target states carry only 28.0 %, indicat- 
ing s evere problems with precision (Kii hner and WhltS 
I1999D . 

There is currently no precise rule to assess the dilemma 
of wanting to target as many states as possible while 
retaining sufficient precision in the description of each 
single one. 

As an example for the excellent performance of this 
method, one may consider the isotropic spin-1 Heisen- 
berg chain, where the single magnon line is shown in 
Fig. El Exact diagonalization, quantum Monte Carlo 
and DMRG are in excellent agreement, with the excep- 
tion of the region g — > 0, where the single- magnon band 
forms only the bottom of a magnon continuum. Here 
Lanczos vector dynamics does not correctly reproduce 
the 2A gap at g = 0, which is much better resolved by 
quantum Monte Carlo. 

The intuition that excitation continua are badly ap- 
proximated by a sum over some O(IO^) effective ex- 
cited states is further corroborated by considering the 
spectral weight function S~^{q = n,uj) [use A — 
in Eq. (|74|l ] for a spin- 1/2 Heisenberg antiferromagnet. 
As shown in Fig. 1151 Lanczos vector dynamics roughly 
catches the right spectral weight, including the 1/w di- 
vergence, as can be seen from the essentially exact cor- 
rection vector curve, but no convergent behavior can be 
observed upon an increase of the number of targeted 
vectors. The very fast Lanczos vector method is thus 
certainly useful to get a quick overview of spectra, but 
not suited to detailed quantitative calculations of excita- 
tion continua, only excitation bands. Nevertheless, this 
method has been applied successfully to the S = h an- 



tiferromagnetic H eisenberg chain (HallberS . ^9^, the 
spin-boson m odel ( Nishi vamal Il99i: ), the Holstein model 
ijZhang et all ri99&). and spin -o rbital chains iii exter- 
nal fields l|Yu and Haasl l2000j) . lOkunishi et~al\ l|200lD 
have used it to ext ract spin chain dispersion relations. 
iGarcia et al\ ^00^ have used continued-fraction tech- 
niques to provide a self-consisten t impurity solver for 
dynamical mean -field calculations l)Georges et all llQQfit 
iMetzner and Vollhardt>.198^ . 



B. Correction vector dynamics 

Even before the advent of DMRG, another way to ob- 
taining more precise spectral functio ns had been pro- 
posed by ISoos and Ramases ha' (* 198 91); it w as firs t ap- 
plied using DMRG bv iRamasesha all l)l997j) and 
iKiihner and Whit3 l)l999|) . After preselection of a fixed 
frequency lu one may introduce a correction vector 



|c(w + i77)) = 



1 



Eq + uj + ir] — H 



i|0), 



which, if known, allows for trivial calculation of the 
Green's function and hence the spectral function at this 
particular frequency: 



GA{uj + iv) = {A\c{uj + iT])). 



(89) 



The correction vector itself is obtained by solving the 
large sparse linear equation system given by 



{Eo+uj + i'n~H)\c{Lu + ir]))=A\0). 



(90) 



To actually solve this nonhermitean equation system, the 
current procedure is to split the correction vector into 
real and imaginary part, to solve the hermitean equation 
for the imaginary part and exploit the relationship to the 
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FIG. 16 Spectral weight ol the S = 1/2 Heisenberg AFM 
from correction vector DMRG. M = 256 states kept. Spec- 
tral weights have been calculated for u-intervals starting from 
various anchoring frequencies for the correction vector. From 
iKiihner and Whita L1999.) . Reprinted with permission. 



real part: 



Re|c(a' + irf}) 



-f]A\Q) (91) 
lm\c{uj + IT])) (92) 



The standard method to solve a large sparse lin- 
ear equation system is the conjugate-gradient method 
llGolub and van Load . Il996() , which effectively generates 
a Krylov space as does the Lanczos algorithm. The 
main implementation work in this method is to pro- 
vide iJ^Imjc). Two remarks are in order. The re- 
duced basis representation of iJ^ is obtained by squar- 
ing the effective Hamiltonian generated by DMRG. This 
approximation is found to work extremely well as long 
as both real and imaginary part of the correction vec- 
tor are included as target vectors: While the real part 
is not needed for the evaluation of spectral functions, 
{Eo + LO - i7)Im|c) Re|c) due to Eq. and target- 
ing Re|c) ensures minimal truncation errors in Hlm\c). 
The fundamental drawback of using a squared Hamil- 
tonian is that for all iterative eigenvalue or equation 
solvers the speed of convergence is determined by the 
matrix condition number which drastically deteriorates 
by the squaring of a matrix. Many schemes are available 
to improve the convergence of conjugate-gradient meth- 
ods, usually based on providing the solution to some re- 
lated, but trivial equation system, such as that formed 
from the diagonal elements of the large sparse matrix 
( Golub and van LoanLll99f]D . 

In the simplest form of the correction vector method, 
the density matrix is formed from targeting four states, 
|0), 1|0), Im|c(cj -f i?7)) an d Re|c(u; -h i?7)). 

As has been shown bv iKiihner and White! l)l999|) . it 
is not necessary to calculate a very dense set of correc- 
tion vectors in w-space to obtain the spectral function 
for an entire frequency interval. Assuming that the fi- 
nite convergence factor -q ensures that an entire range 



of energies of width « is described quite well by the 
correction vector, they have applied the Lanczos vector 
method as detailed in the last section to A|0), using the 
effective Hamiltonian obtained by also targeting the cor- 
rection vector, to obtain the spectral function in some 
interval around their anchor value for oj. Comparing the 
results for some lo obtained by starting from neighbor- 
ing anchoring frequencies allows for excellent convergence 
checks fFig. lTB]) . In fact, they found that best results are 
obtained for a two-correction vector approach where two 
correction vectors are calculated and targeted for two fre- 
quencies ui, L02 = oJi + Aw and the spectral function is 
obtained for the interval [wi,a;2]. This method is, for 
example, able to provide a high precision result for the 
spinon continuum in the S" = 1/2 Heisenberg chain where 
standard Lanczos dynamics fails (Fig. I15|l . 

Reducing the broadening factor 77, it is also possible to 
resolve finite-system peaks in the spectral function, ob- 
taining to some approximation both location and weight 
of the Green's function's poles. 

The correction vector method has been applied to 
determine the nonlinear optical co efficients o f Hub - 
bard chains and de rived models by iPati et aTl (|1999'); 
IKiihner et al have extracted the ac conductivity 

of the B ose-Hubbard mod el with nearest neighbor inter- 
actions. 'Raas et'aR l)2004 h have used it to study the dy- 
namic correlations of single-impurity Anderson models 
and were able to resolve sharp dominant resonances at 
high energies, using optimized algorithms for the matrix 
inversion needed to obtain the correction vector. 



C. Dynamical DMRG 

A further important refinement of DMRG dynamics 
is obtained by a reformulation of the correction vector 
method in terms of a minimization principle, which has 
been called "dynamical DMRG" LJec kclmann, 2002a}. 
While the fundamental approach remains unchanged, the 
large sparse equation system is replaced by a minimiza- 
tion of the functional 

Wa^(w,V)= (93) 
mEo +u-Hf + ry^lV;) + ry(A|V;) + T^i^A). 

At the minimum, the minimizing state is 

IV-min) =Im|c(tj + i?7)). (94) 

Even more importantly, the value of the functional itself 



HOi,r,(w,'0) = -7r??CA(t^ + i?7). 



(95) 



such that for the calculation of the spectral function it is 
not necessary to explicitly use the correction vector. The 
huge advantage is that if the correction vector is known 
to some precision e (which will be essentially identical 
for the equation solver and the minimizer), the value of 
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the functional itself is by general properties of variational 
methods known to the much higher precision e^. Hence 
the DMRG algorithm is essentially implemented as in 
the correction vector method, with the same target vec- 
tors, until they converge under sweeping, but with the 
minimization of WA,riiuJ,tp) replacing the sparse equa- 
tion solver. Results that are obtained for a sequence of 
uji may then be ext e nded t o other uj by suitable interpo- 
lation ijJeckelmannL l2002a(l . also exploiting first deriva- 
tives of spectral functions numerically directly accessible 
at the anchor points. 

For dynamical quantities, there are no strict state- 
ments on convergence. However, convergence for large 
M seems to be monotonic, with Ca typically underesti- 
mated. 

The high-quality numerical data obtained from dy- 
namical DMRG in fact allow for an extrapo lation to the 
thermodynamic limit. As pointed out by iJeckelmannI 
l)2002al) . the double hmit 

C Aiuj) = \im \iia Ca{L;lu + irj), (96) 

where limits may not be interchanged and which is very 
hard to take numerically, may be taken as a single limit 

CaM- lim CAiL;u + irj{L)), (97) 

j)(L)^0 

provided that r]{L) — > as i — > oo is chosen such that 
the finiteness of the system is not visible for the cho- 
sen r]{L) and it thus seems to be in the thermodynamic 
limit from a level spacing perspective. This implies that 
r]{L) > 6uj(L), the maximum level spacing of the finite 
system around energy to. For a typical tight-binding 
Hamiltonian such as the Hubbard model one finds 

viL) > |, (98) 

where c is the bandwidth (and in such Hamiltonians also 
a measure of propagation velocity). The key argument 
is now that if 1](L) is chosen according to that prescrip- 
tion the scaling of numerical results broadened by ri(L) is 
the same as for some thermodynamic limit form known 
or conjectured analytically subject to Lorentzian broad- 
ening using the same rj. From Lorentzian broadening of 
model spectra one can then show that a local (5-peak in an 
otherwide continuous spectrum with weight C'^ scales as 
/iTri[L), and that a power-law divergence [uj — wo)~" 
at a band edge is signalled by a scaling as r? ( L)~°. More 
model spectra are discussed bv I.TeckelmannI l|2002al) . 



Dynamical DMRG has been extensively used to study 
the spectrum and the optical conductivit y of the ex- 



tended Hubbard model ( 
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optical conductivity may be calculated as 



cTi (w) = - — lim ImG J (w -t- i-q) , (99) 
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FIG. 17 Optical conductivity for a Peierls-Hubbard model 
in the [/ » f limit, L = 128, r; = 0.1, n = lo ~ U: dy- 
namical DMRG (solid) vs. broadened exact J-contribution 
(dot-dashed) and unbroadened thermo dynamic hm i t band s 
(dashed). Note log-hnear scale. From Ijeckelmannl i2002ari . 
Reprinted with permission. 

with the current operator J = ^^^-J^ia^ii'^a i'^a i+i ~ 
j^j^Cjj. j). Another possibility is to use 

0-1 (w) = --^ limIm[(w + i7?)Gr,(w + i?7)] , (100) 

with the dipole operator D = —e i{ni — 1). 

In the noninteracting limit U = the optical conduc- 
tivity is nonvanishing in a band 2|A| < to < At with 
square-root edge divergences. Dynamical DMRG repro- 
duces all features of the exact solution including the di- 
vergence quantitatively. It is even possible to extract the 
degree of the singularity. Moving to the strongly inter- 
acting case U ^ t, one expects two continuous bands 
due to the dimerization at 2|A| < \lo — U\ < At and the 
Hubbard resonance at w = U . Fig. 1171 which compares 
to analytical solutions shows both that (up to broad- 
ening at the edges) the bands are reproduced with ex- 
cellent quality and that the ^-singularity is measured 
with extremely high precision. Similarly, very precise 
spectral functions for the Hub bard model away from 
half- filling have been obtained ('Bcn thien et alV l2004|) . 
^ishimoto and Jeckclmann (2004) have shown that dy- 
namical DMRG may also be applied to impurity prob- 
lems such as the flat-band single-impurity Anderson 
model upon suita ble discretization of the band. Recently, 
iNishimoto et al\ l|2004j) have extended this to a precise 
calculation of the Green's function of a single-impurity 
Anderson model with arbitrary band, thereby providing 
a high-quality self-consistent impu rity solver needed for 
dynamical mean-field calcula tions l)Georges et alV Il996t 
iMetzner and Vollhardtlll989fl . 

V. BOSONS AND DMRG 

So far, our discussion of DMRG applications has been 
largely restricted to quantum spin and electronic systems, 
which are characterized by a fixed, usually small number 
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of degrees of freedom per site. As algorithmic perfor- 
mance relies heavily on iVgito small (formally 0{N^^^^), 
in practice rather 0{N^^^^^)), one may wonder whether 
DMRG is applicable to bosonic degrees of freedom with 
iVsitc — oo. Such bosonic degrees of freedom occur for 
example in the Bose-Hubbard model, 



H 



BH 



or fermionic degrees of freedom. While they are be- 
lieved to be reliable to give a generic impression of 
physical phenomena, for more precise studies more ad- 
vanced techniques are nec essary (compare results of 
ICaron and Moukouril l)l996() and of .Bursill (,,1999.') ). 



1), (101) B. Large number of degrees of freedom 



that has come to the forefront of reseach due to the re- 
alization of a tunable quantum phase transition from a 
Mott-insulating to a superfluid pha se in ultracold bosonic 
atomic gases in optical lattices l|Greiner et al\ . l2002|) . 
Another model of interest is the Holstein model, where 
electrons (also spinless fermions or XY spins) couple to 
local (quantum) phonons that react dynamically and are 
not a priori in some adiabatic limit, 

i 

~lY.^b\ + 6,)(n, - 1) + c^^(6l6, + 1/2). 

i i 

It models electrons in a vibrating lattice and opens the 
way to polaron physics. Yet another model is the spin- 
boson model, which models the dissipativc coupling of a 
two-state system to a thermal reservoir given by bosonic 
oscillators: 



HsB — —-^o^ 
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A. Moderate number of degrees of freedom 

The simplest conceptual approach is to arbitrarily 
truncate the local state space to some A'^max for the 
bosonic degrees of freedom, and to check for DMRG con- 
vergence both in M and in A'max- This approach has 
been very successful in the context of the Bose-Hubbard 
model where the onsite Coulomb repulsion /7 sup- 
presses large occupation numbers. It has been used for 
the standard Bose-H u bbard model ijKiihner and Monieii 
Il998t iKiihner el oi.L l200(D . with a random potential 
l|Rapsch et alV Il999j) . in a parab olic potential due to 
a magnetic trap for cold atoms ijKollath et al\ . l2004j) 
and with non-hermitian hopping and a pinning impu- 
rity t o model superconductor flux lines l|Hofstetter et alV 
|2_004') . Typically, allowing for about three to five 
times the average occupation is sufficient as A/'max- 



membranes fNishivama', 


"2002albl 120031) and quantum 


strings ( Nishivama. 




2001) necessitated the introduc- 



tion of a larger number of local vibrational states. 
Other applications that a re more pro blematic have 
been to phonons, both wit h (|CarorL.aiid3Ioukourl Il996t 
iMaurel and Le pctit. 2000:"Mau rel et a/.ll200lfl and with- 
out ijCaron and Moukouri .1997|) coupling to magnetic 



Essentially three approaches have been taken to reduce 
the large, possibly divergent number of states per site to 
a small number manageable by DMRG. 

iBursilll l)l999l) has proposed a so-called jour block ap- 
proach that is particularly suited to periodic boundary 
conditions and is a mixture of Wilson numerical renor- 
malization group and DMRG ideas. Starting from 4 ini- 
tial blocks of size 1 with M states (this may be a rela- 
tively large number of electronic and phononic degrees of 
freedom) forming a ring, one solves for the ground state of 
that Af* state problem; density matrices are then formed 
to project out two blocks, and form a new block of dou- 
ble size with states, which are truncated down to M 
using the density-matrix information. From 4 of these 
blocks, a new 4 block ring is built, leading to a doubling 
of system size at every step. Calculations may be simpli- 
fied beyond the usual time-saving techniques by reducing 
the number of product states to some smaller num- 
ber of states for which the product of their weights in 
the density matrices is in excess of some very small e, 
of the order of 10~^^ to 10~^°. Translation operators 
applied to the resulting ground state of the 4 blocks on 
the ring allow the explicit construction of specific mo- 
mentum eigenstates, i n particu lar for k = and k — n. 
For both an XY-sp in ijBursillL ^99) and isotropic spin 
ijBursill et all Il999j) variety of the Holstein model, this 
method allows to trace out very precisely a Kosterlitz- 
Thouless phase transition from a quasi-long-ranged anti- 
fcrromagnet for small spin-phonon coupling to a dimer- 
ized, gapped phase. As Kosterlitz-Thouless phase tran- 
sitions exhibit an exponential ly slow opening of the gap 
ijOkamoto and Nomural Il992() . the exact localization of 
the transition by analyzing the gap is problematic. It is 
rather convenient to apply the level spectroscopy tech- 
nique l)Nomura and Okamotol Il994|) which locates the 
phase transition by a small-system finite-size extrapola- 
tion of the crossing points g* of the lowest lying excita- 
tion of different symmetries, including k — and k — tt 
states that are easily accessible in the 4 block approach. 
It was found that in the neighborhood of the phase tran- 
sition, roughly 30 phonon states were sufficient to obtain 
converged results. A similar scenario of a KT transition 
was obtained for spinless fermions in the Holstein model 
ijBursill alL Vl99S). 

An approach more in the spirit of DMRG, the so-called 
local s tate reduction, was introduced by IZhang et all 
While it can also be combined with exact di- 
agonalization, I want to formulate it in a DMRG setup. 
Assuming a chain with fermionic and a small number 
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FIG. 18 Optimal vs. bare phonon states: ground state en- 
ergy convergence for a 4-site Holstein model at half-filling vs. 
number of ( o ptimal , bare) phonon states kept. Taken from 
IZhang et al\ |l9983). 



FIG. 19 Standard reorganization of a two-dimensional lattice 
as a zig-zag one- dimensional chain with long-ranged interac- 
tions for DMRG treatment. Typical blocks during a finite- 
system DMRG application are shown. 



of bosonic degrees of freedom on each site, one of the 
sites is chosen as "big site" to which a further number 
of bare bosonic degrees of freedom is added. Within a 
DMRG calculation, a density matrix is formed for the big 
site to truncate the number of degrees of freedom down 
to some fixed number of "optimal" degrees of freedom. 
This procedure is repeated throughout the lattice in the 
finite-system algorithm, sweeping for convergence and for 
adding further bosonic degrees of freedom. The stan- 
dard prediction algorithm makes the calculation quite 
fast. Physical quantities are then extracted within the 
optimal phononic state space. As can be seen from Fig. 
1181 merely keeping two or three optimal states, in which 
high- lying bare states have non-negligible weight, may 
be as efficient as keeping of the order of hundred bare 
states. This approach has allowed to demonstrate the 
strong effect o f quantum lat t ice flu ctuations in trans- 
polyacetylene ijBarford et al I l2002aD . Combined with 
Lanczos vector dynamics, very precise dynamical sus- 
ceptibilities have been extracted for spin-boson models 
("Nishivama, 1999*) . Ext ensions of the method are found 
in .Friedman. (,2000^1 and lFehskd l|2000() . 

iJeckelmann and Whitd l)l998() have devised a further 
approach where 2" bosonic degrees of freedom are em- 
bodied by n fermionic pseudosites: writing the num- 
ber of the bosonic degree of freedom as a binary num- 
ber, the degree of freedom is encoded by empty and full 
fermionic pseudosites as dictated by the binary number. 
All operators on the bosonic degrees of freedom can now 
be translated into (rather complicated) operators on the 
fermionic pseudosites. Finite-system DMRG is then ap- 
plied to the resulting Hamiltonian. They have been able 
to study polaronic self-trapping of electrons in the Hol- 
stein model for up to 128 phonon states and have located 
very precisely the metal- i nsulator transition in this model 
l)Jeckelmann et a]] . ll999j) . 



VI. TWO-DIMENSIONAL QUANTUM SYSTEMS 

Not least since the spectacular discovery of high-Tc su- 
perconductivity related to Cu02 planes, there has been a 
strong focus on two-dimensional quantum systems. Early 
on, it was suggested that the Hubbard model or the t-J 
model away from half-filling might be simple yet power- 
ful enough to capture essential features of high- Tc su- 
perconductivity. The analytical study of these quan- 
tum systems is plagued by similar problems as the one- 
dimensional case, as in many effective field theories d = 2 
is the lower critical dimension and few exact results are 
available. Numerically, exact diagonalization methods 
are even more restricted in two dimensions and quantum 
Monte Carlo is difficult to use for fermionic models away 
from half-filling due to the negative sign problem. Can 
DMRG help? 

The first step in any two-dimensional application of 
DMRG is to identify blocks and sites in order to ap- 
ply the strategies devised in the one-dimensional case. 
Assuming nearest-neighbor interactions on a square lat- 
tice one might organize columns of sites into one super- 
site, such that blocks are built from supersites, sweeping 
through the system horizontally. While this approach 
maintains short-ranged interactions, it must fail for any 
two-dimensional strip of appreciable width L as the num- 
ber of states per supersite grows exponentially as iV^^^g 
and is only useful for narrow ladder systems. 

The standard a p proach ( Noack, White, and S c alapino 
in iL andau et all |l994) . Liang and Panel l)l994j) . IWhitd 
(1996q)), known as multichain approach is to define a 
suitable one-dimensional path over all sites of the two- 
dimensional lattice, such as the configuration shown in 
Fig.lll 

One may now apply standard one-dimensional finite- 
system DMRG at the price of long-ranged interactions of 
range 2L both within and between blocks, as indicated 
in Fig. 119! An inherent difficulty is the preparation of 
blocks and operators for application of the finite-system 
algorithm, as there is no precursor infinite-system DMRG 
run in some sequence of smaller systems possible in this 
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setup. Few compositions of smaller blocks in this scheme 
resem ble the final system at all. While iLiang and Pand 
l|l994l) have simply switched off all non-nearest neigh- 
bor interactions in the mapped one-dimensional system 
and applied standard infinite-system DMRG in order to 
switch on all interactions in finite-system runs, one can 
also grow blocks in a standard infinite-system DMRG 
run, where some very short exactly solvable system are 
used as environment such that the entire superblock re- 
mains treatable. This and similar warmup procedures 
will generate starting states for finite-system DMRG 
far from anything physically realistic. As finite-system 
DMRG provides only sequential local updates to that 
state, there is no guarantee that the system does not get 
trapped in some local minimum state. That such trap- 
pings do exist is seen from frequent observations that 
seemingly converged systems may, after many further 
sweeps, suddenly experience major "ground state" en- 
ergy drops into some new (possibly) ground state. 

White and coworkers have followed the approach to ex- 
ploit local trappings by theorizing ahead of using DMRG 
on possible ground state types using other analytical or 
numerical techniques and to force these types of states 
onto the system by the application of suitable local 
magnetic fields and chemical potentials. These exter- 
nal fields are then switched off and the convergence be- 
havior of the competing proposed states under finite- 
system DMRG observed. This procedure generates a 
set of states each of which corresponds to a local en- 
ergy minimum. The associated ph ysical properties may 
be mea sured and compared (see IWhite and Scalapind 
l)l998a|) and [White and ScalaoinO (200(f) for a discus- 
sion). In the multi chain approach , the two-dimensional 
Heisenberg model fWhitd. Il996bl). the two-dimensional 
t-J-model fChernvchev et al. '. '2003t iKamnf et all l200lt 
^hitc and Scalaoino, 1998a b, 2000) with particular em- 
phasis on stripe formation, and the 6- leg ladder Hubbard 
model ijWhite and Scalapind f2003|) have been studied. 

Am ong competing setups llFarnell|. [2003: HenehuJ 
1999t Idu Croo dc Jona h and va njjgeuwenl. ll99S 
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Moukouri and Caron. ■2003 V 'Xian g et al\ l|200H) have 



set up a two-dimensional algorithm that uses a true 
DMRG calculation all along and builds L x L systems 
from previously generated {L — 1) x {L — 1) systems 
while keeping the lattice structure intact. It can be 
applied to all lattices that can be arranged to be square 
with suitable interactions: e.g., a triangular lattice is 
a square lattice with additional next-nearest neighbor 
interactions along one diagonal in each square. 

A one-dimensional path is organized as shown in Fig. 
1201 where the pair of free site may be zipped through the 
square along the path using the standard finite-system 
algorithm, yielding arbitrary block sizes. In particular, 
one may obtain triangular upper or lower blocks as shown 
in Fig.|^ Combining these blocks from a (L — 1) x (L — 
1) system and adding two free sites at the corners, one 
arrives at a L x L system. Here, the upper left free site 
can be zipped to be a neighbor of the lower right free site. 




system 



FIG. 20 Diag onal re organization of a two-dimensional lattice 
as used bv lXiang et ah (200 li) . Typical blocks during a finite- 
system DMRG application are shown. 
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FIG. 21 Composition of blocks from a {L—l)x{L—l) system 
and 2 free sites into a, L x L system as used by|Xiane et gl| 
ll200lf) . The fat solid line indicates the one-dimensional path 
through the lattice. The lattice subsets surrounded by the 
broken lines are blocks A and B; the last added sites are in- 
dicated by broken circles. Note that blocks overlap for the 
smaller lattice and are "pulled apart" for the big lattice. The 
open and full circles stand for the new sites added. 



as it sits next to active block ends (i.e. the ends where 
new sites are added). The pair of free sites can now be 
zipped through the system to yield the desired triangular 
blocks for the step L ^ L + 1. 

Both for two-dimensional square and triangular S = 
1/2 Heisenberg models this approach systematically 
yields lower energies than the multichain approach with 
the exception of very small systems. Even fo r a relatively 
modes t number of states kept (M = 300) IXiang et all 
l)200lj) report thermodynamic limit extrapolations in 



good agreement with large-scale quantum Monte Carlo 
simulations: For the square lattice, they find Eoo = 
-0.3346 versus a QMC result Er^ = -0.334719(3) 
l|Sandvikl Il997|) . The potential of this approach seems 
far from exploited at the moment. 

So far, I have been tacitly assuming that interactions 
are of the same order along both spatial directions. Ex- 
perimentally, a relatively frequent situation is that one- 
dimensional quantum systems are weakly coupled along 
a second axis. At high enough temperatures this interac- 
tion will washed out, but at sufficiently low temperatures 
there will be a crossover from one- to two-dimensional 
behavior. Such systems may be studied by attempt- 
ing a precise one-dimensional description and introduc- 
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ing the weak interchain interacti on on the mean field 
or som e ot her perturbative level. iMoukouri and CaronI 
l|200.'^ and lMoukouTil l|20n4l) have used DMRG in a simi- 
lar spirit by solving a one-dimensional system using stan- 
dard DMRG with M states, and determining the M' 
lowest-lying states for the superblock. Chain lengths are 
chosen such that the lowest lying excitation of the finite 
chain is down to the order of the interchain coupling. 
These states are taken to be the states of a "site" , which 
are combined to a new effective Hamiltonian, that then 
is once again treated using DMRG. Results for weakly 
coupled spin chains are in good agreement with quan- 
tum Monte Carlo results, however M' is severely limited 
to currently several tens. 

A severe limitation on two-dimensional DMRG is pro- 
vided by the exponential growth of M with L for a pre- 
selected truncated weight or ground state precision (Sec. 
IIILBjl . For frustrated and fermionic systems beyond the 
very small exact diagonalization range (6 x 6 for S* = 1/2 
spins) DMRG may yet be the method of choice as quan- 
tum Monte Carlo suffers from the negative sign prob- 
lem and even medium-sized fermionic systems of size, 
say, 12 X 12 sites would be most useful; in models with 
non-abelian symmetries, the implementation of the asso- 
ciated good quantum numbers has been shown to reduce 
drasticall y the truncation error for a given number of 
state s (.McCuUoch et all l200ll: iMcCulloch and Gulacsl 

1200(1 l2nniLl2on2D. 

Very recently, IVerstraete and Cirad l)2004|) have pro- 
posed a new approach to two-dimensional quantum sys- 
tems combining a generali zation of their matrix p rod- 
uct formulation of DMRG (Verstrae te et ad l2004bll and 
imaginary-time evolution (,Verstraete''er'fflIf l2004a|^ . dis- 
cussed in Sec. IIII.AI and Sec. IIX.CI One-dimensional 
matrix-product state are formed from matrix products of 
M X M matrices (^[tT])^^, with A/-dimensional auxiliary 
state spaces on the bond to the left and right of each site. 
The two-dimensional generalization for a square lattice is 
given by the contraction of tensors {A[a])af3'yS with M- 
dimensional auxiliary state spaces on the four adjacent 
bonds. Finite-temperature and ground states are calcu- 
lated by imaginar y-time evolution o f a tota lly mixed state 
along the lines of IVerstraete et As tensorial 

contractions lead (unlike in the case of the ansatz matri- 
ces A[a] appearing in one dimension) to a proliferation of 
indices on resulting tensors, a suita ble truncat ion scheme 
is needed as described bv .Verstraete and Cirac 1,2004) . 

An appealing feature of this approach is that the en- 
tropy of entanglement for a cut through a system of size 
L X L will scale, for fixed M, linearly in system size as 
it should generically. The errors in energy seem to de- 
crease exponentially in M. M is currently very small (up 
to 5), as the algorithm scales badly in M; however, as M 
variational parameters are available on each bond, even 
a small M corresponds to large variational spaces. At 
the time of writing, it is too early to assess the potential 
of this new approach; however, current results and its 
conceptual clarity make it seem very promising to me. 



VII. BEYOND REAL SPACE LATTICES 

In this section, I shall consider three groups of DMRG 
applications that seem to have very little in common at 
first sight: the study of translationally invariant low- 
dimensional models in momentum space, high-precision 
quantum chemistry calculations for small molecules, and 
studies of small grains and nuclei. However, from a 
DMRG point of view, all share fundamental properties. 
Let me first discuss momentum-space DMRG, move on to 
quantum chemistry and finish by considering small grain 
and nuclear physics. 

A. Momentum-space DMRG 

Real-space DMRG precision dramatically deteriorates 
when applied to long-ranged interactions or hoppings. 
Moreover, momentum is not accessible as good quantum 
number in real-space DMRG. Momentum-space DMRG, 
on the other hand, makes momentum a good quantum 
number, works naturally with periodic boundary con- 
ditions, and allows trivial manipulation of the single- 
particle dispersion relation. Momentum-space DMRG 
has already yielded highly satisfying dispersion relations 
and momentum distributions in excellent agreement with 
analytical results (Nishimoto et ai, 2002a). 

Consider for definiteness a Hubbard model with local 
interaction U, but long-ranged hopping tij, 

Hlr = J2 "^T"*i (104) 

ija i 

in arbitrary dimensions; site i is localized at and = 
t{Yi ~ Vj). With L lattice sites, Fourier transformations 

yield the band structure 

6(k)=^e-"'-i(r) (106) 

r 

and the Hamiltonian 

^ = ^'^^^yl.}S'^<yM^Y; ^ 4ki 4k2'^ik3^Tki+k2 -ka- 
le, (T ki,k2.k3 

(107) 

The reciprocal lattice points k are now taken to be 
"sites". Note that L is not the linear size of the lat- 
tice, but is taken to be the total lattice size (because of 
the typical DMRG mapping to one dimension). 

The key idea (Xiang, 1996) is to arrange these "sites" 
into a one-dimensional chain with long-ranged interac- 
tions w hich is trea ted by the real-spac e finit e-system 
DMRG ijLeeeza et a i. 2003a; Nishimoto 'e^ a^J . |2p02a). 
In addition to the conventional good quantum numbers, 
states will also be classified by total momentum, which 
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allows a further substantial segmentation of Hilbert space 
(by a factor of order L), hence for the same number of 
states kept, momentum-space DMRG is much faster than 
real-space DMRG. 

To obtain an efficient implementation, several issues 
have to be addressed. 

(i) In momentum space there is a huge proliferation of 
interaction terms of O(L^), that have to be generated, 
stored and applied to wave functions efficiently. 

(ii) For an application of the finite-system DMRG al- 
gorithm, we need to provide blocks and operators living 
on these blocks for all block sizes. In real space, they are 
naturally generated during the infinite-system stage of 
the algorithm. In momentum space, another "warm-up" 
procedure must be found. 

(iii) On a one-dimensional lattice with short-ranged in- 
teractions it is natural to group sites as they are arranged 
in real space. In momentum space with long-ranged in- 
teractions, there is no obvious "site" sequence, even in 
the one-dimensional case. Does the arrangement affect 
convergence properties, is there an optimal arrangement? 

Let us discuss these points in the following. 

(i) DMRG operator representation for general two-body 
interaction Hamiltonians. In principle we should like to 
use DMRG to treat the generic Hamiltonian 

-fftwo-body = J2 "^IS + ^VklcUyk<^l ' (10^) 
ij ijkl 

where encodes kinetic energy and one-body potentials 
and Vijki a generic two-body (Coulomb) interaction; for 
simplicity, we assume spin contained in the orbital in- 
dices. In the worst-case scenario, all Vijki are different, 
although symmetries and model properties may yield de- 
cisive simplifications. In momentum space, for example, 
Vijki 7^ for one I only once ijk is fixed due to momen- 
tum conservation. 

As for operators living on the same block 

{m\clcj\rh) ^ ^(m|c||m')(m'|c^.|m), (109) 

rn' 

all operator pairings have to be stored separately. Leav- 
ing aside the simpler case where one or two of the four 
operators live on the single sites in DMRG, and assuming 
that they are all either in block S or E, this suggests that 
for L "sites" of the order of 0{L'^) operators have to be 
stored. As their form changes for each of the L possi- 
ble block sizes, the total memory consumption would be 
0{L^AP) on disk for all blocks and 0{L^AP) in RAM 
for the current block. At the same time, for each of the 
L steps of a sweep, calculation time would be of order 
0{L'^M^), or 0{L^M^) for the entire calculation. 

Memory consumption as well as the associated calcu- 
lation time can however be reduced drastically (Xiang, 
119961) . Let us consider the three possible operator distri- 
butions on blocks S and E. 

(a) Four operators in one block (4/0): Terms 
Vijkic\c^jCy.Ci are absorbed into a single block Hamilto- 
nian operator during block growth. Assuming i, j, k 



are in the previous block and site I is added to form the 
current block, a representation of c|cjcj, in the previous 

block basis allows to form Viju x clcjC^. x C; in the block 
plus site product basis, which is then transformed into 
the basis of the current block and added into the single 
block Hamiltonian operator. For L blocks, 0{L^) repre- 
sentations each of c\c^-Ci^ are necessary. These, in turn 
can be compounded into complementary operators 

Oi ^Y.^^,Mc\c]c^, (110) 

ijk 

SO that 

Y.'^^3Mc\c]c,c^^Y.0lCl. (111) 

ijkl I 

The complementary operators can be constructed as dis- 
cussed in Sec. III.GI assuming the knowledge of two- 
operator terms c|cj. For L blocks, 0{L'^) of those exist, 
leading to memory consumption 0{L^M^) on disk and 
0{L^AP) in RAM. 

(b) Three operators in a block (3/1): One applies the 
strategy of (|110|l and Hlll|) . with Oi and C; acting on 
different blocks. 

(c) Two operators in a block (2/2): Again, the com- 
plementary operator technique can be applied, with the 
modification that each complementary operator living on 
block S has now two matching operators in E. A further 
class of complementary operators 

Oki=Y.^^3kiclcl (112) 

allows the simplification 

J2 "^vkicUhk^^i ■ (113) 

ijkl kl 

Memory consumption for the second type of complemen- 
tary operator is 0{L^AP) on disk and 0{L'^M'^) in RAM. 
Taking all operator combinations together, global mem- 
ory consumption is to leading order 0{L^M^) on disk 
and 0{L'^M'^) in RAM, which is a reduction by com- 
pared to the naive estimate. In momentum space, due 
to momentum conservation, memory consumption is re- 
duced by another factor of L to 0{L?KP) on disk and 
0{LM^) in RAM. 

Using the complementary operator technique, calcula- 
tion times are dominated by the (2/2) terms. In analogy 
to the construction of (4/0) terms "on the ffy" by gen- 
erating the new terms of the sum, transforming them 
and adding them into the operator, the Oki can be con- 
structed at a computational expense of 0{L^M^) for gen- 
erating the L new terms to be added to each of the 
complementary operators with Af ^ matrix elements each 
and O(L^M^) for transforming the operators into 
the current basis. Using the multiplication technique 
of Sec. III. II multiplying the Hamiltonian to the state 



33 



vector costs O(L^M^) time at each step or 0{L^M^) 
per sweep. Global calculation time per sweep is thus 
0{L^M^) + 0{L^AP), a reduction by for the domi- 
nant first term (typically, M ^ L for the relevant DMRG 
applications). 

(ii) Setting up finite-system blocks. The stan- 
dard practise currently adopted in momentum space 
l|Nishimoto et a/.l.l2002atlXiand.[T99 6|) is to use standard 
Wilsonian renormalization llWilsonL 119751 119831) : for the 
chosen sequence of momenta, blocks are grown linearly 
(as in infinite-system DMRG), but diagonalized without 
superblock embedding and the lowest-energy states re- 
tained, with the important modification that it should be 
ensured that at least one or two states are retained for 
each relevant Hilbert space sector, i.e. all sectors within 
a certain spread about the average momentum, parti- 
cle number and magnetization. Otherwise DMRG may 
miss the right ground state, as certain sectors that make 
an important contribution to the ground state may not 
be constructed during the finite-system sweep. Left and 
right block sectors must be chosen such that all of them 
find a partner in the other block to combine to the de- 
sired total momentum, particle number and magnetiza- 
tion. Moreover, it is of advantage to choose correspond- 
ing states in different sectors such that obvious symme- 
tries are not broken, e.g. ^ —S^ if the final state is 
known to be in the S^^f = sector; it has also been found 
that choosing initial block sizes too small may result in 
much slower converg ence or trapping in wrong minima 
ijLegeza et aZ.L l2003a|) . 

(iii) Determin i ng the order of momentum "sites". 
iNishimoto et all l)2002aj) report that for short-ranged 
Hubbard models ordering the levels according to in- 
creasing distance to the Fermi energy of the non- 
interacting limit, |e(k) — ep\, thus grouping levels that 
have strong scattering together, works best, whereas 
for longer-ranged Hubbard models a simple ordering ac- 
cording to e (k) works bes t. Similar results have been 
reported by iLeeeza et all (2003a) , who could demon- 
strat e trapping in wrong minim a for inadequate order- 
ings; iLeeeza and Solvoml l)2003() have provided a quan- 
tum information based method to avoid such orderings. 

INishimoto et Ml ^l2002a^ have carried out extensive 
convergence studies for Hubbard models in ID and 2D 
with various hoppings. Convergence seems to be dom- 
inated by the ratio of interaction to bandwidth, U/W 
rather than some U/t. The momentum-space DMRG 
algorithm becomes exact in the U/W — > limit; conver- 
gence deteriorates drastically with U/W, but is accept- 
able for U/W < 1, which is U < 8t for the 2D Hub- 
bard model with nearest-neighbor hopping. Generally it 
seems that for momentum-space DMRG convergence de- 
pends only weakly on the range of hopping as opposed 
to real-space DMRG. Momentum-space DMRG is worst 
at half-filling and somewhat deteriorates with dimension 
(if calculations for the same physical scale U/W are com- 
pared). In two dimensions, for moderate U, momentum- 
space DMRG is more efficient than real space DMRG. 



Extr apolation schemes. As c an be seen from the re- 
sults of lNishimoto et all l)2002af) , even for values of M as 
large as several 1000, convergence is not achieved. Due 
to the complex growth scheme, it cannot be taken for 
granted even after many sweeps that the environmental 
error has been eliminated; as in two-dimensional real- 
space DMRG, sudden drops in energy are observed. In 
their application, Nishimoto et al. (2002a,'l report that for 
fixed M a fit formula in l/M, 



E 



fit 



M 



a-i 
AP 
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works very well and improve results by over an order of 
magnitude even for final M = 2000. Building on the pro- 
portionality between truncated we ight and energy error 
for eliminated environmental error, ILeeeza et aHl2003aj) 
have proposed to vary M to maintain a fixed truncated 
weight ("dynamical block selection scheme"). They find 
that in that approach 



Ejep) ~ E, 



LaCn 



(115) 



is very well satisfied for a « 1, where Cp is the truncated 
weight. A further advantage of that approach is that 
the number of states needed for a certain precision varies 
widely in momentum-space DMRG such that for many 
DMRG steps important savings in calculation time can 
be realized. 

All in all, momentum-space DMRG seems a promis- 
ing alternative to real-space DMRG, in particular for 
longer-ranged interactions (for short-ranged interactions 
in one dimension, real-space DMRG remains the method 
of choice). However, momentum-space DMRG presents 
additional complications (optimal ordering of levels, effi- 
cient encoding) that may not have been solved optimally 
yet. 



B. Quantum chemistry 

A field where DMRG will make increasingly important 
contributions in the next few years is quantum chemistry. 
While the first qu antum-chem i cal DM RG calculations on 
cy clic polyene by iFano et all l)l99^ and polyacetylene 
by iBendazzoli et all l)l999fl were still very much in the 
spirit of extended Hubbard models, more recent work has 
moved on to calculations in generic bases with arbitrary 
interactions. Let me situate this latter kind of DMRG 
applications in the general context of quantum chemistry. 

Within the framework of the Born-Oppenheimer ap- 
proximation, i.e. fixed nuclear positions, two major ways 
of determining the electronic properties of molecules are 
given by Hartree-Fock (HF) and post-HF calculations 
and density functional theory (DFT). DFT is computa- 
tionally rather inexpensive and well-suited to quick and 
quite reliable studies of medium-sized molecules, but not 
overly precise in particular for small molecules. Here, 
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Hartree-Fock calculations, incorporating Fermi statistics 
exactly and electron-electron interactions on a mean-field 
level provide good (initial) results and a starting point 
for further refinement, so-called post-HF calculations. 
Quantum chemistry DMRG is one such post-HF calcula- 
tion. 

For the HF and post-HF calculations, one starts from 
a first-quantized Schrodinger equation (in atomic units) 
for the full molecule with N electrons; second quantiza- 
tion is achieved by introducing some suitably chosen basis 
set In this way, a general two-body interaction 

Hamiltonian as in Eq. pUHII can be derived. 

A HF calculation now solves this Hamiltonian at mean- 
field level, and it can be reexpressed in terms of the HF 
orbitals. A closed-shell singlet ground state is then sim- 
ply given by 

|HF) =cji4i...c|^/2cj^/2|0). (116) 

Within standard quantum chemistry post-HF calcu- 
lations, one improves now on this ground state by a 
plethora of methods. Simple approaches consist in diag- 
onalizing in the subspace of one- or two-particle excited 
states; others choose a certain subset of electrons and di- 
agonalize the Hamiltonian fully in one subset of orbitals 
for those (complete active space). Taking all electrons 
(excitations) and all orbitals into account, one obtains, 
within the originally chosen basis set, the exact many- 
body solution, known as the fuU-CI (configuration inter- 
action) solution. While the approximate schemes typi- 
cally scale as to iV^ in the number of orbitals, fuU- 
CI solutions demand exponential effort and are available 
only for some very small molecules. 

From the DMRG point of view, determining the 
ground state of the full second-quantized Hamiltonian is 
formally equivalent to an application of the momentum- 
space DMRG, with some additional complications. As 
momentum conservation does not hold for Vijki, there 
are 0{N) more operators. Moreover, the basis set must 
be chosen carefully from a quantum chemist's point of 
view. As in momentum space DMRG, a good sequence 
of the strongly different orbitals has to be established for 
DMRG treatment such that convergence is optimized. 
Also, an initial setup of blocks must be provided, includ- 
ing operators. 

Orbital ordering in quantum chemistry turns out to 
be crucial to the performance of the method; a badly 
chosen ordering may lead to much slower convergence 
or even trapping in some higher energy local mini- 
mum. The simplest information available, the HF en- 
erg y and occupation o f the o r bitals, has been exp loited 
bv iWhite and Martini lll999l). iDaul et all ()2000a|l . and 
iMitrushe nkov et al\ ll200ll |2003|) to enforce various or- 
dcrings of the orbitals. 

iGhan and Head-GordonI (2002) have proposed to use 
the symmetric reverse Cuthill-McKcc (CM K) reorder - 
ing CCuthill and McKcc, 1969; Liu and Shcr mail Il975|l . 
where orbitals are reordered such that Ty matrix ele- 
ments above a certain threshold are as band-diagonal as 



possible. This reduces effective interaction range, leading 
to faster con vergence in M and red uced sweep number, as 
confirmed bv lLegeza et al\ l)2003b|) . who optimize the po- 
sition of the Fock term + occupied (^^ifc/cj - "^Vikjk) 
and report reductions in M by a third or so. 

Another approach is given by exploiting the short 
range of chemical interactions which in the basis of typ- 
ically delocalized HF orbitals leads to a large effective 
range of interactions detrimental to DMRG. This may 
be reduced by changing to more localized orbitals. Or- 
bital locali z ation t echniques have first been studied by 
iDaul et a/] (l2000a ) who a pplied a localization procedure 
(jPipe^'andMezev , Il989|) to the unoccupied orbitals to 
avoid increases in the starting configuration energy, but 
did not find remarkable improvement over orderings of 
the unmodified HF orbitals. 

Quantum inform a tion t echniques have been used by 
iLeeeza and Solvoml l)2003l) both for momentum space 
and quantum chemistry DMRG to devise optimal order- 
ings. Studying various ad hoc orderings, they find fastest 
and most stable convergence under sweeping for those 
orderings that maximize entanglement entropy [Eq. H22|l ] 
between orbitals and the rest of the chain for orbitals at 
the chain center; orderings are hence chosen such that 
strongly entangled orbitals are brought to the chain cen- 
ter (often those closest to the Fermi surface). This or- 
dering can be constructed iteratively starting from a first 
guess as obtained by the CMK ordering. 

Once the orbital ordering has been selected, current 
applications follow a variety of "warm-up" procedures to 
build all sizes of blocks and to ensure that block states 
contain as many states as possible relevant for the fi- 
nal result. This may be done by augmenting the block 
being built by a toy environment block with some low 
energy states that provides an environment compatible 
with the block, i.e. leading to the desired total quantum 
numbers JCh an and Head- GordonLl200l . doing infinite- 
system DMRG on left and right blocks while targeting 
also excited states (Mitrushen kov et ai, 2 001) , and defin- 
ing a minimum number of states to be kept even if the 
non-zero eigenvalue s of the density m atrix do not provide 
enough states (Leg eza et all l2003a|) . Block state choice 
can also be guided by enta nglement entropy calculations 
l|Legeza and Solvoml 1200,^ . 

Various molecules have by now been studied using 
DMRG, both for equilibrium and out-of equilibrium con- 
figurations, ground states and excitations. H2O has been 
serving as a benchmark in many quantum chemistry cal- 
culations and has been studied within various finite one- 
particle bases, more recently in the so-called DZP basis 
llBauschlicher and Tavloiill986() with 8 electrons in 25 or- 
bitals (and 2 "frozen" electrons in the Is oxygen orbital) 
and the larger TZ2P basis fWidm ark et ai, 1990) of 10 
electrons in 41 orbitals. While the larger basis will gen- 
erally yield lower energies, one can compare how close 
various approximations come to an exact diagonaliza- 
tion (full configuration interaction (CI) - if at all pos- 
sible) solution within one particular basis set. In the 
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smaller DZP basis, the exact ground state energy is at 
-76,256634H (Hartree). While the Hartree-Fock solu- 
tion is several hundred mH above, and the single and 
doubles configuration interaction solution about ten mH 
too high, and variou s coupled cluster approximations 
l|Bartlett et all Il99(i) are reaching chemical accuracy of 
1 mH with errors between 4.1 mH and 0.2 mH, DMRG 
results of various groups l |C han and He ad-Gordon. 2002; 
iDaul all l2000at IWhite and Martini 11999^1 are better 
than 0.2 mH at M ~ 4 00, rea ching an error of about 
0.004 mH at M ^ 900 ijChan a nd Hcad-Gordo^, l2002t 
iLeeeza aILl2003a^. Moving to the larger TZ2P basis 
llChan and Head- GordonL 12001 ■ where there is no exact 
solution, accuracies of all methods are ranking similarly, 
with DMRG outperforming the best coupled cluster re- 
sults with an error of 0.019 mH at M ~ 3000 and reaching 
an error of 0.005 mH at AI ^ 5000, with the extrapolated 
result being — 76.314715H. The comparison is even more 
favorable if the nuclei are arranged away from their equi- 
librium configuration. 

Excellent performance is also observed in 
the extra ction of dissocia ti on e nergy curves 
for N2 llMitrushenkov aZ.I . l200l|) and water 
l)Mitrushenkov et all l200^ 

Low-lying excited states have been studied for HHeH 
I Paul et ad l2000al). CH9 ijLeeeza et all l2003aj) and LiF 



i Legeza et oi.l . l2003b() : the latter case has provided a test 
system to study how efficiently methods map out the 
ionic-neutral crossover of this alcali-halogenide for in- 
creasing bond length. With relatively modest numerical 
effort (M not in excess of 500 even at the computation- 
ally most expensive avoided crossing) relative errors of 
10~^, 10"'', and 10~^ for the ground state energy, first 
excited state energy and the dipole moment compared 
to the full CI solution l)Bauschlicher and LanghoM Il988|) 
have been obtained. The dipole moment precision in- 
creases by orders of magnitude away from the avoided 
crossing. 

To assess the further potential of quantum chemistry 
DMRG for larger molecules and orbital bases, it should 
be kept in mind that the generic scaling of the algo- 
rithm - which does compare extremely favorably with 
competing methods that scale as to N'^ - is modified 
by several factors. On the one hand, for a given desired 
error in energy, various authors report, as for conven- 
tional DMRG, that 6E ^ Ep, the truncated weight, which 
scales (Sec. IHI!b|| as ep — exp(— fcln^M). k, however, 
shrinks with size in a model-dependent fashion, the most 
favorable case being that orbitals have essentially chain- 
like interactions with a few "nearest neighbor" orbitals 
only. It may be more realistic to think about molecular 
orbitals as arranged on a quasi one-dimensional strip of 
some width L related to the number of locally interact- 
ing orbitals; in that case, for standard strip geometries, 
the constant has been found to scale as L~^. So iV* scal- 
ing is in my view overly optimistic. On the other hand, 
on larger molecules orbital localization techniques may 
be more powerful than on the small molecules studied 



so far, such that orbital interactions become much more 
sparse and scaling may actually improve. 

One other possible way of improving the scaling prop- 
erties of quantum chemistry DMRG might be the canon- 
ical diagonalization approach proposed bv I Whit'3 lf2002) . 
which attempts to transform away numerically by a se- 
quence of canonical basis transformations as many of the 
non-diagonal matrix elements of the second-quantized 
Hamiltonian H of Eq. (|108|l as possible such that en- 
tire orbitals can be removed from H, resulting in a new, 
smaller quantum chemistry problem in a different ba- 
sis set, which may then be attacked by some DMRG 
technique such as those outlined above. The removed 
orbitals, which are no longer identical with the original 
ones, are typically strongly overlapping with the energet- 
ically very high-lying and low-lying orbitals of the origi- 
nal problem that are almost filled or empty. Of course, 
these transformations cannot be carried out for the ex- 
ponen tially large number of Hilbert space states. IWhitel 
l)2002|) moves to the HF basis and carries out a particle- 
hole transformation; the canonical transformations are 
then carried out in the much smaller space of states not 
annihilated by one of the Vijki and formed by the mini- 
mum of creation operators from the HF vacuum (ground 
state). For example, the term Vijkidldjdkdi, where the 
dl,d.i are particle-hole transformed fermionic operators, 
implies the consideration of the "left" state dj|0) with HF 
energy Ei and the "right" state d]„|0) with HF en- 
ergy Er- In this smaller state space, a sequence of canoni- 
cal basis transformations may be impleme nted as a differ- 
ential equation as originally proposed bvlWeenerj (Il994) 
as flow equation method and iGlazek and WilsonI |)1994|) 
as similarity renormalization: with A some antihermitian 
operator, a formal time-dependence is introduced to the 
unitary transformation H{t) — exp(fyl)iJ (0) exp(— i^), 
where H{Q) is the original Hamiltonian of Eq. H108|) ex- 
pressed in the HF basis. The corresponding differential 
equation 



dH{t) 
dt 



[A.H{t)] 



(117) 



is modified by making A time-dependent itself. One now 
expands 



and 



i7(t) = ^a„(t)/io 



A{t) = y^5aaa(t)faa. 



(118) 



(119) 



where each ha is the product of creation and annihilation 
operators and each Sq some constant yet to be chosen 
under the constraint of the anti-hermiticity of A. To 
avoid operator algebra. White introduces 



[K.hp] =^Blph^, 



(120) 
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where additional contributions to the commutator that 
cannot be expressed by one of the operator terms in Eq. 
H118|l are suppressed; this eUminates three-body inter- 
actions generically generated by the transformation and 
not present in the original Hamiltonian. Then Eq. (|117|1 . 
with A time-dependent, reduces to a set of differential 
equations 

^=J2Blp^^^a{t)a,{t), (121) 

a/3 

that can now be integrated numerically up to some time 
t. It can now be shown that the goal of diminishing aa{t) 
for all non-diagonal contributions to H can be achieved 
efficiently by setting 

s^^iEi-Ery\ (122) 

where Ei and Er are the HE energies of the left and right 
HF orbitals in ha. Stopping after some finite time t, one 
can now remove orbitals (e.g. with the lowest occupancy, 
i.e. almost full or empty before the particle-hole transfor- 
mation) and use DMRG for the full many- body problem 
in the reduced orbital space. IWhitd l)2002|) finds this ap- 
proach to yield very good results as compared to fuU-CI 
calculations for a stretched water molecule, but as of now 
the full potential of the method remains unexplored. 

To conclude, in my opinion, DMRG has been estab- 
lished as a competitive, essentially fuU-CI quality method 
in quantum chemistry which should be taken seriously. 
However, efficient geometrical optimization techniques to 
determine lowest-energy molecular geometries are still 
lacking. 

C. DMRG for small grains and nuclear physics 

Let us reconsider the generic Hamiltonian of Eq. 
(|108|l . where energy levels above and below the Eermi 
energy are ordered in ascending fashion and where a 
simple model interaction is chosen. This Hamiltonian 
can then be treated in the following so-called particle- 
hole reformulation of infinite-system real-space DMRG 
(j^ukclskv an d Sierra..l999. 2000): Starting from two ini- 
tial blocks from a small number of states, one ("hole 
block") from the lowest unoccupied one-particle levels 
above the Eermi energy and the other ( "particle block" ) 
from the highest occupied one-particle levels below, we 
iteratively add further states to particle and hole block, 
moving away from the Eermi surface, determine the 
ground state, and in the usual DMRG procedure cal- 
culate density matrices and the reduced basis transfor- 
mations by tracing out particle and hole states respec- 
tively. This approach may work if there is a clear physical 
reason that states far away from the Fermi surface will 
hardly influence low energy physics and if interactions 
show "translational invariance" in energy space; other- 
wise the infinite-system algorithm should fail. The ad- 
vantage of such simple models is that much larger sys- 
tems can be treated than in quantum chemistry. 



A model Hamiltonian which happens to combine both 
features is the reduced BCS Hamiltonian 

Hbcs = Y^^ej - i^)cl^Cja "^"^Yl ''lAi^n^'n ' (12^) 

jo- ij 

where DMRG has been able to definitely set- 
tle longstanding questions llDukelskv and Pittei l200lt 
'Duke lskv and SierraL Il999l l2000^ on the nature of the 
breakdown of BCS superconductivity in small grains ex- 
pected when the level spacing d in the fin ite system 
becom es of the order of the BCS gap A l)Andersonl 
Il959j) . DMRG has conclusively shown that there is a 
smooth crossover between a normal and a superconduct- 
ing regime in the pairing order parameter and other 
quantities, in contradiction to analytical approaches in- 
dicating a sharp crossover. This approach to DMRG has 
been extended to the study of the Josephson effect be- 
tween superconducting nan ograins with discrete energy 
levels (iGobert et all l20n4bD and to the observation and 
calculation of well-defined quasi particle excitations in 
small interacting disor dered systems with h igh dimen- 
sionless conductance g l|Gobert et 0^11200431) . 

Particle-hole DMRG has a lso been applie d suc- 
cessfully in nuclea r physics JPimitrova et all l2002t 
iDukelskv and Pittei 120011: iDukelskv et all l2002() in the 
framework of the nuclear shell model, where a nucleus 
is modeled by core orbitals completely filled with neu- 
trons and protons ("doubly magic core") and valence 
orbitals partially filled with nuclcons. The core is con- 
sidered inert and, starting from some Hartree-Eock level 
orbital configuration for the "valence" nucleons, a two- 
body Hamiltonian for these nucleons is solved using the 
particle-hole method. This is conceptually very similar 
to the post-HE approaches of quantum chemistry, with 
model interactions for the nucleons as they are not so pre- 
cisely known, such as pairing interactions, quadrupolar 
interactions, and the like. 

This approach has been very successful for nuclcons in 
very large angular momentum shells interacting through 
a pairing and a quadrupola r force in an oblate nucleus 
touk elskv and PittelL 12001)) . with up to 20 particles of 
angular momentum j = 55/2, obtaining energies con- 
verged up to 10~^ for M = 60; for 40 nucleons of j = 
99/2, 4 digit precision was possible; in t his case, 38 % of 
energ y was contained in the correlations ifDukelskv et all 
|2002|) . For more realistic calculations of the nucleus 
^^Mg, with 4 neutrons and 4 protons outside the dou- 
bly magic ^^O core, convergence in M was so slow that 
almost the complete Hilbert space had to be considered 
for good convergence ijDimitrova et all l2002|) . The fun- 
damental drawback of the particle-hole approach is that 
it does not allow for an easy implementation of the finite- 
system algorithm (levels far away from the Eermi surface 
hardly couple to the system, giving little relevant infor- 
mation via the density matrix) an d that angular mom en- 
tum conservation is not exploited. iPittel et all l)2003l) are 
currently aiming at implementing an algorithm for nuclei 
using angular momentum, circumventing these difhcul- 
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ties (see also iDukelskv and Pittej l)200 



VIII. TRANSFER MATRIX DMRG 

Conventional DMRG is essentially restricted to 
T = calculations, with some computationally ex- 
pensive fo rays into the very-low temperature re gimes 
possible llBatista e.t nlV Il998t iHallberg et lUX Il999t 
iMoukouri and Ca ronlll996|) . Decisive progress was made 
by Nishino tl995|} who realized that the DMRG idea 
could also be used for the decimation of transfer ma- 
trices, leading to the name of transfer-matrix renormal- 
ization group (TMRG). This opened the way to DMRG 
studies of classical statistical mechanics at finite temper- 
ature for systems on two-dimensional L x oo strips. If 
one applies the generic mapping of d-dimensional quan- 
tum systems at finite temperature to {d+ l)-dimensional 
classical systems, TMRG also permits study of thermo- 
dynamic properties of one-dimensional quantum systems 
at finite temperature. 



A. Classical 2D transfer matrix DMRG: TMRG 

Consider the textbook transfer matrix method for a 
one-dimensional classical system with Hamiltonian H 
and A'^gite states Icr,) on each site i, e.g. the Ising model. 
Assuming nearest-neighbor interactions, one performs a 
local decomposition of the Hamiltonian H = hi, where 
each hi contains the interaction between sites i and i + 
Taking the partition function of a system of length N 
with periodic boundary conditions, 



=Tre-'5«-Tre-Sti/3''. 



(124) 



and inserting identities |(Ti)(cri|, one obtains for a 
translationally invariant system ~ Tr T^, where the 
TVsitc X iVsitc transfer matrix T reads 



(a,\T\a 



i+l/ 



-/3hi 



From this one deduces for N ^ oo the free energy per site 
/ = —kBT\nXi[T], where Ai is the largest eigenvalue of 
T, from which all thermodynamic properties follow. 

Moving now to a two-dimensional geometry, on a strip 
of dimension L x N with a short-ranged Hamiltonian, 
the transfer matrix method may be applied in the limit 
N ^ oo for some finite width L ^ N, treating a row of 
L sites as one big site; transfer matrix size then grows as 
iVj'jg, strongly limiting this approach in complete analogy 
to the exact diagonalization of quantum Hamiltonians. 
Hence (Fig. I22|l the first dimension, with finite system 
length L, is attacked by DMRG applied to the trans- 
fer instead of the Hamiltonian matrix, to keep transfer 
matrix size finite and sufficiently small by decimation, 
while retaining the most important information. Results 
are extrapolated to infinite size. The second dimension. 
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FIG. 22 Strategy for the DMRG treatment of two- 
dimensional classical systems. 



with infinite system length, is accounted for by the one- 
dimensional transfer matrix method. 

To prove that this concept works, i.e. that an optimal 
decimation principle can be set up for transfer matrices 
in analogy to the case made for Hamiltonians, consider a 
short-ranged classical Hamiltonian atr>OonaLxAf 
strip, where N oo. We now define an unnormalized 
density-matrix pu — e~^^ (in reality an evolution oper- 
ator) by 



(125) 



where a and cr label states of the L sites on top and 
bottom of the strip. T^^^ is the band transfer matrix of 
Fig. [52 The partition function Z = Tr /5„ is then given 
as 



E 



(126) 



where Ai > A2 > . . . are the eigenvalues of T^^^f; the 
largest being positive and non-degenerate for positive 
Boltzmann weights, partition function and unnormal- 
ized density-matrix simplify in the thermodynamic limit 
TV cx) to 



Af 



and 



/5„ = Af |Ai)(Ai|, 



(127) 



(128) 



where |Ai) is the normalized eigenvector of the largest 
transfer matrix eigenvalue. 

Consider now Fig. 1231 The unnormalized density- 
matrix has exactly the same form, up to the prefactor, 
of the pure state projector of standard DMRG, with |Ai) 
assuming the role of the target state there. Tracing out 
the right "environment" half, the left "system" unnor- 
malized reduced density-matrix is given by 



N 



L/2 



PuS = Tr^Af |Ai)(Ai| = Af ^ w^\w^){wa.\, (129) 
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FIG. 23 Pictorial representation of the partition function, 
the (unnormalized) density matrix and the reduced density- 
matrix for a two-dimensional classical system on a strip ge- 
ometry. The black stripe represents a band tr ansfer matrix. 
Adapted from Nishino in lPeschel et al 



where Wa are the positive eigenvalues to the normalized 
eigenvectors \wa) of the reduced "system" density ma- 
trix; Wa — 1. Completing the partial trace, one finds 



a=l 



(130) 



such that indeed the best approximation to Z is obtained 
by retaining in a reduced basis the eigenvectors to the 
largest eigenvalues Wa as in conventional DMRG. 

Let us explain the DMRG transfer-matrix renormaliza- 
tion (jNishino. ,.1995,,) in more detail. Because the transfer 
matrix factorizes into local transfer matrices, some mi- 
nor modifications of the standard DMRG approach are 
convenient. Assume as given a M x M transfer matrix 
of a system of length L, 



{ma 



|r( 



(131) 



where the sites at the active (growth) end are kept ex- 
plicitly as opposed to standard DMRG, and |m), \rh) are 
block states fFisf. I24|) . Considering now the superblock 
of length 2L, the transfer matrix reads 



(132) 



Lanczos diagonalization yields the eigenvector of T^^^^ 
to the maximum eigenvalue Ai. Again in some analogy 
to conventional DMRG, the fact that the transfer ma- 
trix is a product can be used to speed up calculations by 
decomposing |(/)) — T^^^^\^) into three successive multi- 
plications, 10) = r(^)[r(2)(rW|V>))]. 

The obtained eigenvector now yields the reduced den- 
sity matrices and the reduced basis transformations. The 
last step is to build the Af A'sito x MNsHc transfer matrix 
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FIG. 24 Transfer matrix DMRG step. 
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for a system of length L + 1 (Fis. I25|l: 



(133) 



J2 Hnr) {nT\T^^'> |nf ) {Ta\T^^'>\fa) (nf |m) . 

nhrf 

The TMRG procedure is repeated up to the desired 
final length, the free energy then obtained from Ai; this 
allows the calculation of all thermodynamic quantities 
through numerical differentiation of the free energy. To 
avoid numerically unstable second derivatives for spe- 
cific heat and susceptibility, it is convenient to consider 
the first order derivative of the average energy extracted 
from expectation values such as (cTiCTj) = (Ai jCTiCTj |Ai) 

obtained by replacements e^'"'' — > die^^^^ in the above 
calculations; at the same time, correlation lengths may 
be extracted both from these expectation values or from 
the two leading transfer matrix eigenvalues. 



e = -l/lnRc(A2/Ai). 



(134) 



To refine results, a finite-system calculation can be set 
up as in conventional DMRG. A main technical differ- 
ence between conventional DMRG and classical transfer 
matrix DMRG is the absence of good quantum numbers, 
which simplifies the algorithm, but is hard on computa- 
tional resources. 
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A large body of works has emerged that actually 
strongly relies on the finite strip-width provided by the 
TMRG, studying competing bulk and surface eS^ects both 
in two-dimensional Ising models on a L x oo strip. The 
confined Ising model has been used to model two coex- 
isting phases in the presence of bulk and surface fields as 
well as gravity, in order to model wetting and coexistence 
phenomena and the competition of bulk and surface ef- 
fects in finite geometries. In the case of opposing surface 
fields favoring phase coexistence at zero bulk field up to 
some temperature that (unint uitively) goes to the wet - 
ting temperature for L ^ oo ijParrv and Evanst ll99Cl|) . 
it could be shown that gravity along the finite width di- 
rection suppresses fluctuations such that it restores the 
bulk critic al temperature as limiting tem peratu re for co- 
existence l|Carlon and Drzewinski . 19971 Il998|) . These 
studies were extended to the compet ition of surface and 
bulk fields if Przewihski et al\. Il998|) . Further studies 
elucidated finite-size corrections to the Kelvin equation 
at complete wetting for parallel surface and bulk fields 
(ICarlon e t al\ |l99 8tl. the nature of coexisting excited 
state s Jp^ewiAskil l2000l). the sc a ling o f cumulant ra- 
tios l|Drzewinski and Woitkiewic3 . l200(]|) . and an anal- 
ysis of the confined Ising model at and near criticality 
for cam peting surface and bulk fields l|Maciolek e.t ol\ . 
I2001alb^ . In the case of the Potts model with Q > A, 
where there is a first-order phase transition in the bulk 
l|Baxteil[T98l. but a seco nd-order phase transition at the 
surface ( Lipowskvl Il982() , TMRG permitted to extract 
surface critical exponents that demonstrate that the uni- 
versality class of the sur face transition is Q-independent 
l|lgloi and Carlonill999D . 

TMRG has also been used to study criti cal properties 
of trul y two-dimensional classi c a,l syst e ms ifCarlon ei~ai. , 
1999al: [Honda and Horiguchi', Il997l : ISato and Sasakil 
20001: [Tsushima et al, 1997), such as the spin-3/2 



Ising model on a square lattice ((Tsushim^^^ 



^997*), the Ising model with line defects (jChung et 
[2000), the three-state chiral clock model as exam- 
ple of a commensurate- incommensurate phase transition 
llSato and Sasakil I2fl0n|) . the 19-vertex model in two di- 
mensions as a discrete version of the continuous XY 
mode l to study t he Kosterlitz-Thouless phase transi- 
tion l|Honda and Horiguchi. 1997.) , and the random ex- 
change coupling classical Q-state Potts model which was 
demonstrated to have critical properties independent 
of Q (Carlon et al, 1999a). Lav and Rudnick (200j) 
have used the TMRG to study a continuous Ginzburg- 
Landau field-theory formulation of the two-dimensional 
Ising model by retaining the largest eigenvalue eigenfunc- 
tions of bond transfer matrices as state space. 

iGendiar and Surdal l)2000D have extended the TMRG 
to periodic boundary conditions to obtain, at the expense 
of the expected lower DMRG precision, much better ther- 
modynamic limit scaling behavior for L — > 00. Critical 
temperature, thermal and magnetic critical exponent are 
all extracted with modest computational effort at a stag- 
gering 6 to 7 digit precision compared to the Onsager 
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FIG. 26 Corner transfer matrix setup for the calculation 
of reduced density matrices and partition functions in two- 
dimensional classical models. 
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B. Corner transfer-matrix DMRG 



Following iBaxted figil, one may conclude that the 
essential aspect of the reduced density-matrix is that it 
represents a half-cut in the setup of Fig. 1231 whose spatial 
extension is to be taken to the thermodynamic limit. The 
same setup can be obtained considering the geometry of 
Fig. [^ where four corner transfer-matrices ( CTM) are 
defined as 



(<t'|C( 
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Cr int 
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{ijki) 



a., a. 



|r( 



2) 
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where the product runs over all site plaquettes and the 
sum over all internal site configurations, i.e. the states 
on the sites linking to neighboring CTMs are kept fixed. 
The corner site state is invariant under application of the 
CTM. The unnormalized reduced density-matrix is then 
given by 



E 

cr(i)(T(3)cr(2) 



(*t(5)|c(^)|^(4))(^(4)|cW|^(3)) 

(*t(3)|C(^)|^(2))(^(2)|cW|^(1)) 



(136) 



with the center site fixed. Diagonalizing C^^-*, and ob- 
taining eigenvalues and eigenvectors |Ai), the reduced 
density-matrix Pu{L) then reads 



UL)^Y.^t\\^){X^ 



(137) 



and the partition function, obtained by tracing out the 
reduced density matrix, is 



(138) 



Following the ar gu ment for classical TMRG, 
iNishino and OkunishJ l|l996[l have introduced the 
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FIG. 27 Corner transfer matrix growth using band and pla- 
quette transfer matrices. Solid states are summed over. 



CTMRG or corner transfer-matrix renormalization 
group. A sequence of increasingly large corner transfer 
matrices is built (Fig. I27|l as 



(139) 



'ma) 



where previous reduced basis transformations arc sup- 
posed to have taken place for C'^^-' and T'^^. Simi- 
larly, r(^+i) is built as in classical TMRG. Diagonalizing 
(7(^+1) yields eigenvalues and associated eigenvectors, 
the M most important of which define the reduced ba- 
sis transformation. This reduced basis transformation is 
then carried out on T'^^^^^ and (7*^^+^^ completing one 
CTMRG step. The crucial advantage of this algorithm 
is that it completely avoids the diagonalization of some 
large sparse matrix. 

CTMRG allows to calculate similar quantities as 
TMRG, all thermodynamic variables as some derivative 
of the free energy or using local e xpectation values. For 
the tw o-dimensional Ising model iNishino and Okunishil 
l)l997j) obtained Ising critical exponents at 4 digit pre- 
cision, using systems of up to L = 20000. For the 
(5 = 5 two-dimensional Potts model, which has a very 
weak first-order transition they could determine the 
latent heat £ « 0.027 compared to an exact C — 
0.0265 llBaxteil [l982), which is hard to see in Monte 
Carlo simulations due to metastability problems. Other 
applications have conside r ed th e spin-3/2 Ising model 
l)Tsushima and Horiguchil Il998^ and a vertex model 
with 7 vertex con figurations, mod eling an order-disorder 
transition (^Takasaki et all l200l|) . More recent exten- 
sions study self-avoiding walk models in two dimensions 
l|Foster and Pinettesl l2003alb^ . 

One may also generalize th e CTMRG to three- 
dime nsional classical models (jNishino and Okunishl 
Il998|) . While the implementation is cumbersome, the 
idea is simple: the semi-infinite cut line in the plane lead- 
ing to 4 corner matrices gives way to a semi-infinite cut 
plane in some volume leading to 8 corner tensors two 



of which have an "open" side. Growing these corner 
tensors involves concurrent growing of corner matrices, 
band transfer matrices, and plaquette transfer matrices. 
Numerical results, however, indicate clear practical limi- 
tations of the concept. 



C. Quantum transfer matrix DMRG in ID 

It is now straightforward - at least in principle - to 
extend classical TMRG to the calculation of the thermo- 
dynamics of one-dimensional quantum systems due to the 
general relationship between d-dimensional quantum and 
{d + l)-dimension al classical problem s. This approach 
was first used by Bursill et al. fl996') and full y devel- 
oped by IWa ng an d XLaag. (1997J as well as by iShibatal 
l)l997j) . To appreciate the following, it is best to visu- 
alize quantum TMRG as an extension of the q uantum 
transfer-matrix method (Bctsuvaku, 1981 Il985|) in the 
same way conventional DMRG extends exact diagonal- 
ization. 

Consider the mapping of one-dimensional quantum to 
two-dimensional classical systems, as used in t he quan- 
tum transfer matrix and quantum Monte Carlo ijSuzukiL 
Il976j) methods: Assume a Hamiltonian that contains 
nearest-neighbor interactions only and is invariant un- 
der translations by an even number of sites. We now 
introduce the well-known checkerboard decomposition 
lISuzTikJ. IT97I H = H1 + H2. Hi = EfJi h2^-l and 
H2 = ^'ii- is tti^ local Hamiltonian linking sites 

i and i + l; neighboring local Hamiltonians will in general 
not commute, hence [Hi,H2] 7^ 0. However, all terms in 
Hi or H2 commute. N is the system size in real space 
(which we w ill take to infin ity). 

Following 'Trottei^ fl959*), we consider a sequence of 
approximate partition functions 



Zl ■■= Tr (e 



with Trotter number L. Then 



-IJHi/L -13H2/L^L 



Z = lim Zl 



(140) 



(141) 



holds: quantum effects are suppressed as neglected com- 
mutators between local Hamiltonians scale as (/3/L)^. 

One now expands Hi and H2 in the exponentials of 
Zr m Eq. l(nnjl into the sums of local Hamiltonians and 
inserts 2L + 1 times the identity decomposition 



N 



(142) 



sandwiching each exponential. Introducing an addi- 
tional label j for the identity decompositions, we find 
(2L -I- 1) X "sites", where each site carries two labels, 
corresponding to two dimensions, the real space coordi- 
nate (subscript) i and the Trotter space (imaginary time) 
coordinate (superscript) j. Thinking of both coordinates 
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FIG. 28 Checkerboard decomposition: active vs. inactive pla- 
quettes of the two-dimensional effective classical model. 



as spatial, one obtains a two dimensional lattice of dimen- 
sion (2L -I- 1) X A^, for which the approximate partition 
function Zi^ reads: 



N/2 

Trn 



2j + l 2J4-1 
'2i-l "21 



-Ph2 



I 2?" 21 



"2i+l / 



Figure [2H1 shows that as in a checkerboard only every sec- 
ond plaquette of the two-dimensional lattice is active (i.e. 
black). To evaluate the trace is taken over all local 
states of the (2L-I-1) x N sites, while noting that the trace 
in Eq. H14U|) enforces periodic boundary conditions along 



the imaginary time direction, \a]) = \a. 



2L+1\ 



Nothing 



specific needs to be assumed about boundary conditions 
along the real space direction. Note that the orientation 
of the transfer matrix has changed compared to Fig. 1221 
Working the way backward from the partition func- 
tion Zl, we can now identify transfer and density ma- 
trices. Introducing a local transfer-matrix as — 
exp{— (Jhk / L) , and exploiting the assumed restricted 
translational invariance, the global transfer matrices 

L 

]-J(^2j + 1^2, + l|^^|^2,^2j^ 

L 

J = l 

can be defined (see Fig. 1281) , representing odd and even 
bonds on the chain because of the alternating checker- 
board decomposition. The local transfer matrices are 
linking the states of two sites at different Trotter times 
and are evaluated using the Trotter-time independent 
eigenbasis of the local Hamiltonian, h\i) = Ei\i): 



-PE./L 



(143) 



Summing over all internal degrees of freedom 
obtain matrix elements 



we 



a \I1I2\T ...T 



2L + 1\ 



(144) 



Using the global transfer 
density-matrix reads 



matrices, the unnormalizcd 



and 



PuL = [T1T2] 



Zl = Tr [T1T2] 



N/2 



N/2 



(145) 



(146) 



somewhat modified from the classical TMRG expression. 

As Zl — Af ^^ + A^^^ + . . ., where Ai are the eigenvalues 
of 7i72, the largest eigenvalue Ai of 7i72 dominates in 
the thermodynamic limit iV — > 00. The density matrix 
simplifies to 



PuL = Af/^|V'fl)(V'L|, 



(147) 



where (V'lI and \iPr) are the left and right eigenvectors to 
eigenvalue Ai. Due to the checkerboard decomposition, 
the transfer matrix 7172 is non-symmetric, such that left 
and right eigenvectors are not identical. Normalization 
to (V'lIV'-r) = 1 is assumed. The free energy per site is 
given as 



ifcsTlnAi. 



(148) 



If Zl (i.e. the eigenvalue Ai) is calculated exactly, com- 
putational effort scales exponentially with L as in exact 
diagonalization. As the errors scale as (/3/I/)^, reliable 
calculations are restricted to small /3, and the interesting 
low-temperature limit is inaccessible. Instead, one can 
adopt the classical TMRG on the checkerboard transfer 
matrix to access large L, hence low temperatures. 

Conceptually, the main steps of the algorithm are un- 
changed from the classical case and the argument for 
the optimality of the decimation procedure is unchanged. 
Explicit construction and decimation formulae are how- 
ever slightly changed from the classical case and much 
more complicated notatio nally because of the che cker- 
board decomposition : see IWang and Xian and 
IShibatal Ill997ll2003^ . 

At the start, the transfer matrix T^^^ of two plaquettes 
is given by (Fig. 



J ^(CTlCTa, (71(73,(72^3) = 

E ^'^^ ^^2 I n I (7lV|) (a| all f2 I (7 1 (7 1 ) . 



(149) 



The addition of plaquettes follows a zig-zag pattern. 
Let me just discuss the case where the transfer matrix 
grows to comprise an even number of plaquettes. This 
number be L. Then the growth formula reads fFig. I29() : 



42 



L+1 
1 

<P Q 

r1 



L+1 
2 



^2 



3 3 
"2 

Q © 



r1 



© — e 

4 < 



-9 



© — © 

°2 °: 



O © 

1 

°2 



^L+1 

after transformation 



(a) 



(b) 



FIG. 29 Two plaquette transfer matrix T'^-* and transfer ma- 
trix growth. Solid circles are summed over. 



I distinguish between inner and outer states, the outer 
states being those of the first and last row of the transfer 
matrix. For the inner states, the compound notation mi 
and 7713 indicates that they may have been subject to a 
reduced basis transformation beforehand. If the number 
of states miaf and msai^ is in excess of the number 
of states to be kept, a reduced basis transformation is 
carried out, micrf rhi and msCTg — > 7713, using the 
reduced basis transformation as obtained in the previous 
superblock diagonalization. 

The superblock transfer matrix of 2L plaquettes is now 
formed from the product of two such T(^+^\ summing 
over states = o\^^^ and ct^+^ (Fig. OOI). This su- 
perblock transfer matrix is diagonalized using some large 
sparse eigensolver to get the maximum eigenvalue Ai and 
the right and left eigenvectors (V'lI and IV'fl), which may 
be chosen mutually biorthonormal; (V'lIV'-r) — 1- In view 
of the possible need for a reduced basis transformation 
in the next step, non-symmetric reduced density matrices 
are formed. 



/5 = Tr |V'i?,)(^z 



(150) 



where trace is over all states in columns 1 and 3 except 



m , 771 , CT]^ and CTg as they are the states subject to 
a subsequent reduced basis transformation. Row and col- 
umn matrices of the M left and right eigenvectors with 
highest eigenvalue weight yield the transformation ma- 
trices for left and right basis states, i.e. in columns 1 and 
3. This construction implies that basis vectors (jn\ \ and 
\'m-i) will also be biorthonormal. The procedure is re- 
peated until the system has reached the desired final size 
of 2Lijiax plaquettes. 

Several technical issues have to be discussed for quan- 
tum TMRG. An important reduction of the computa- 
tional load comes from the existence of modified good 
quantum numbers derived from conservation laws of 
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FIG. 30 Quantum transfer matrix DMRG step. Solid circles 
are summed over; note the periodic boundary conditions. 



the underlying Hamiltonian. Their existence was al- 
ready observed in the quantu m transfer matrix method 
ijNomura and YamadaL Il99l[l : The action of a local 
magnetization-conserving Hamiltonian in Trotter di- 
rection between imaginary times j and i + 1 obeys 

\s{^^Y + \si^^Y = \s{Y + \Si\\ or, - \s{Y - 

-{S^^^Y + This generalizes to Y.j{~^y\S{Y = 

— ^^(— 1)^[S'2]^. Thus the staggered magnetization, 
summed in Trotter direction, is a conserved quantity, i.e. 
constant along the real-space axis: 



2L 

E 



= V(_l)»+^[5|]^=cst. 



(151) 



In fermionic systems, particle number conservation im- 



plies that P = (— l)'"*"-^?!^ is an add itional good 

quantum number ijShibata and TsunetsueuL Il999a|) . As 
in T = DMRG the "good quantum numbers" are con- 
served by the DMRG decimation process. Moreover, the 
eigenstate with the maximum eigenvalue is always in the 
subspace Q = (or P = 0, respectively): Thermody- 
namic quantities are independent of boundary conditions 
in high-temperature phases, which are found for all finite 
temperatures in one dimension. For open boundary con- 
ditions in real space, the lack of an update at the open 
ends at every second Trotter time step implies the sub- 
traction of equal quantities when forming Q, leading to 
Q = 0. 

An important technical complication occurs because 
the transfer matrices to be diagonalized are asymmet- 
ric due to the checkerboard decomposition. We there- 
fore have to distinguish between left and right eigenvec- 
tors (t/^lI and It/'a). Nonhermitian diagonalization rou- 
tines lack the stabilizing variational principle of the her- 
mitian case and have to be monitored extremely care- 
fully. The simplest choice, very slow but stable, is the 
power method applied to the transfer matrix and its 
transpose. Faster algorithms that have been used to 
obtain stable results are the Arnoldi and unsy mmet- 
ric Lanczos methods l|Golub and van Loanl Il996|) . com- 
bined with rebiorthogonalization. Once the eigenvectors 
have been obtained, the unsymmetric density-matrix [Eq. 
(|150|l ] has to be fully diagonalized. Here, two numerical 
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problems typically occur, appearance of complex conju- 
gate eigenvalue pairs with spurious imaginary parts and 
loss of biorthonormality. The former is dealt with by 
discarding the spurious imaginary parts and by taking 
the real and imaginary part of one of the eigenvectors 
as two real-valued eigenvectors. The latter, less fre- 
quent one, is dealt with by iterative reorthogonalization 
jAmmon et al, 1999). 

Both left and right eigenvectors and {iPr) are 

needed to construct the density matrix as in Eq. (|150|l . 
Having both eigenvectors, a symmetric choice (with 



/Ssymm = Tr^ (^fl | + (152) 



for the density matrix is also conceivable, if one argues 
that this amounts to trying to optimally represent both 
the left and right eigenvectors at the same time. This 
choice is much more stable numerically, as no complex 
eigenvalue problem appears and diagonalization routines 
are much more stable, and it somehow also accounts for 
the infor mation carried by both eigen vectors. A detailed 
study by iNishino and Shibatal l)l999l) clearly favors the 
asymmetric choice. They showed that if the asymmetric 
density- matrix has real eigenvalues, it is for M retained 
states as precise as the symmetric choice for 2M states; 
if they are complex, there is no such advantage. 

For the extraction of physical quantities it is now of 
advantage not to fix the inverse temperature f3 and vary 
the Trotter number L, but to replace /3/L by a fixed ini- 
tial inverse temperature /So ^ 1 in Eq. (|140|l . Quantum 
TMRG growth then reaches for 2L plaquettes all inverse 
temperatures /3o, 2/3o, . . . , LPol at each step we obtain the 
free energy per site /(T) and thus all thermodynamic 
quantities, such as the internal energy u, magnetization 
TO, specific heat at constant volume Cy, and magnetic 
susceptibility x by numerical differentiation. Temper- 
atures less than a hundredth of the main energy scale 
of the Hamiltonian can be reached. Convergence must 
be checked both in M ^ cxd and in /3o ^ 0; the latter 
convergence is in /3g. For a fixed number of plaquettes, 
finite-system DMRG runs can be carried out. 

Severe numerical problems may occur for the second 
order derivatives c„ and x- A direct calculation of u and 
TO, which are expectation values, reduces the number of 
differentiations to one. 

To do this, consider the internal energy (u) = 
Z^^Ti [hie^^^]. In the Trotter decomposition Ti 
changes to: 



vectors of Ai) 



(153) 



and therefore {{tpL\ and \iPr) are the left (right) eigen- 



Tr [e-f"] 



Al 



(154) 



A direct calculation of Cy from energy fluctuations is less 
accurate. The susceptibility can be obtained from calcu- 
lating the magnetization m at two infinitesimally differ- 
ent fields. 

In the beginning, applications of the quantum TMRG 
focused on the thermodynamics of simple Heisenberg spin 
chain s l)ShibataL Il997l: IWang and Xiand. Il997l: IXiantd. 
Il998|) : but modified (e.g., frustrated and dimerized) 



spin ch ains are also easily acc essibl e (U ohnston et al 
' 200(i iKaradamodou et all l200lt iKliimpcr et al 
19991 l200d: iMaeshiina and Okunishj I200C : 



Maisineer and SchoUwofi^i ITflfl^ Ishibata and Ueda , 
2D01; Wang et al.. 20Qft) ; for frustrated chains, effective 
sites of two original sites are formed to ensure nearest- 
neighbor interactions and the applicability. Quantum 
TMRG studies of ferrimagnetic spin chains have revealed 
an interesting combination of thermo dynamics typical 
of ferromagnets and antiferromagnets l | K olezhu k ei^ all 
ll999HMaisinger et ad ll998HYamamoto et all\l99^ . In 
fact, the results are remarkably similar to thos e obtained 
for spin ladders ((Hagiwara et al.. . .2000; .Wang and YuL 
120001) . 

Electronic degrees of freedom have been studied 
for the t-J mode l bv l Ammon et all l)l999|) , and 
ISirker and Kliimperl ll2002albt) and for the Kondo lattic e 
ijShibata et all Il998tlshibata and Tsunetsugul Il999bfl . 
Lately, quantum TMRG has also been applied to 
the algorithmically very s imilar spin-orbit chain by 
ISirker and KhaliuUi J ('2003'). 

Rommcr and Egg crt ( 1999j studied one localized im- 
purity linking two spin chains in the thermodynamic 
limit; there has also been strong inte rest in the case of 
multiple mobile impurities jAmm on~a nd Imada, 2000a ^ 

Dynamic properties at finite temperature are among 
the most frequently available experimental data. Quan- 
tum TMRG offers, in close analogy to quantum Monte 
Carlo calculations, a way to calculate these properties, 
albeit with far less than the usual precision. It has b een 
applied to anisotropic spin chains ("Naef et al., "1999'), to 
the one-dimensional Kondo chain (^Mutou et al. . 1998bl 
ll999l:IShibata and Tsunetsugul Hgbgaj) a nd has been used 
to ca lculate nuclear relaxation rates ijNaef and Wangl 
I2OOOI) . One starts by calculating imaginary-time correla- 
tions G(t) by inserting two operators at different imagi- 
nary times and time distance r into the transfer matrix 
analogous to Eq. p53|l . The spectral function A{uj) can 
now be extracted from the well-known relationship to 
G(r), 



1 r 



duj K{t,uj)A{uj), 



(155) 
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where the kernel K{t,uj) is 

K{t, uj) = e-™ + 6-'^'^+™ (156) 

and the extension to negative frequencies obtained 
through A{—u;) = e~^'^A{uj). This decompo sition is not 
uniqu e; authors also use related expressions ijNaef et al\ . 
^93), but the essential difficulty in using these expres- 
sions is invariant: If we introduce a finite number of dis- 
crete imaginary times ti and frequencies Wj , and a set 
Kij = K{Ti,ujj), Gi = G{Ti), Aj = A{ujj) the spectral 
function A{lo) is in principle obtained from inverting (by 
carrying out a singular value decomposition of Kij) 

G,=Y,K^,A,, (157) 

j 

which, as has been known to practitioners of quantum 
Monte Carlo for a long time, is a numerically unstable 
procedure, because is very badly conditioned. The 
ratio of its largest eigenvalues to its smallest eigenvalues 
is normally so big that it is not e ncoded properly in com- 
puter number representations. iNaef et all 1)1999(1 have 
argued, as is done in quantum Monte Carlo, that given 
the algorithm-based imprecisions in the imaginary-time 
data, one may merely try to find that spectral function 
that is (in a probabilistic sense) most compatible with 
the raw data. This maximum entropy procedure leads to 
finding that set {A{ujj)} that maximizes 

aS[A]~]pcM, (158) 

where S[A\ is the relative entropy 

S[A] = ^(^, - m, - A, ln(A,/m,))(l + e"'^-) (159) 
j 

and [A\ accommodates the noisiness of the raw data, 

i 3 

In the relative entropy, rrij is some spectral function that 
embodies previous knowledge relative to which informa- 
tion of A{ujj) is measured. It may simply be assumed flat. 
The factor 1 -I- e~^^ ensures that contributions for w < 
are considered correctly, af measures the estimated error 
of the raw data and is in TMRG of the order of 10~^. 
It is found by looking for that order of magnitude of af 
where results change least under varying it. a is deter- 
mined self-consistently such that of all solutions that can 
be generated for various a the o ne most probable in the 
light of the input data is chosen l|Jarrell and Gubernatid . 
11996). 

Another way of calculating A{lu) is given by Fourier 
transforming the raw imaginary time data to Matsubara 
frequencies on the imaginary axis, w„ — {2n + 1)/Tr for 
fermions and cj„ = 2n/7r for bosons, using = i(3/L and 

G(ic.„)-^Ee-"-'G(TO. (161) 

1=0 



G(iujn) is written in some Pade approximation and then 
analytically continued from ia;„ to infinitesimally above 
the real frequency axis, lu + ii], j] = 0^ . 

The precision of quantum TMRG dynamics is far lower 
than that of T = dynamics. Essential spectral features 
are captured, but results are not very reliable quanti- 
tatively. A serious alternative is currently emerging in 
time-dependent DMRG methods at finite temperature 
fSec lIX-mi . but their potential in this field has not been 
demonstrated yet. 



IX. SYSTEMS OUT OF EQUILIBRIUM 

The study of strongly correlated electronic systems, 
from which DMRG originated, is dominated by attempts 
to understand equilibrium properties of these systems. 
Hence, the theoretical framework underlying all previ- 
ous considerations has been that of equilibrium statisti- 
cal mechanics, which has been extremely well established 
for a long time. Much less understanding and in partic- 
ular no unifying framework is available so far for non- 
equilibrium physical systems, the main difficulty being 
the absence of a canonical (Gibbs-Boltzmann) ensemble. 
Physical questions to be studied range from transport 
through quantum dots with strong voltage bias in the 
leads far from linear response, to reaction-diffusion pro- 
cesses in chemistry, to ion transport in biological sys- 
tems, to traffic flow on multilane systems. Quite generi- 
cally these situations may be described by time-evolution 
rules that, if cast in s ome operator la nguage, lead to non- 
hermitian operators ('Glaubeil. ll963j) . 

The application of DMRG to such problems in on e 
spatial dimension was pioneered by iHieidal l|l998f) . 
who studied, using a transfer-matrix approach, the 
discrete-time asymmetric exclusion process (ASEP), 
a biased hopping of hard-core particles on a chain 
with injection and removal of particles at both ends 
l|Derrida et all Il993tl . Excellent agreement with exact 
solutions for particular parameter sets (which correspond 
to matrix-pr oduct ansatzes) and with other numerics 
was found. iKaulke and Pesclieil l|l998|) studied the q- 
symmetric Heisenberg model out of equilibrium. 



A. Transition matrices 

In the field of non-equilibrium statistical mechanics 
there has been special interest in those one-dimensional 
systems with transitions from active to absorbing steady 
states under a change of parameters. Active steady 
states show internal particle dynamics, whereas absorb- 
ing steady states are truly frozen (e.g. a vacuum state). 
These transitions are characterized by universal critical 
exponents, quite analogously to equilibriu m phase tran- 
sitions' see Hinrichsen ( 2000) for a review. ICarlon et al\ 
l|l999bj) have shown that DMRG is able to provide reli- 
able estimates of critical exponents for non-equilibrium 
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phase transitions in one-dimensional reaction-diffusion 
models. To this end they applied DMRG to a (non- 
hermitian) master equation of a pseudo-Schrodinger 
form. This approach has been at the basis of most later 
DMRG work on out-of-equilibrium phenomena. Consid- 
ering a chain of length L, in which each site is either 
occupied by a particle (A) or empty (0), the time evolu- 
tion of the system is given in terms of microscopic rules 
involving only neighboring sites. Typical rules of evolu- 
tion are site-to-site hopping (diffusion) yl0 or the 
annihilation of neighboring particles A A 00, all with 
some fixed rates. Non-equilibrium phase transitions in 
the steady state may arise for competing reactions, such 
as A A 00 and v40, 0A A A simultaneously present. 
The last model is referred to as the contact process. Its 
steady-state transition between steady states of finite and 
vanishing density in the thermodynamic limit is in the 
universality class of directed percolation. 

Once the reaction rates are given, the stochastic evolu- 
tion follows from a master equation, which can be written 
as 

= -H\Pit)), (162) 

where \P{t)) is the state vector containing the probabili- 
ties of all configurations, as opposed to the true quantum 
case, where the coefficients are probability amplitudes. 
The elements of the "Hamiltonian" H (which is rather a 
transition matrix) are given by 

(a\H\T) = — w(r a) ; a ^ t 

{a\H\a) = J2w{a^T), (163) 

where \a), \t) are the state vectors of two particle con- 
figurations and w(t — *■ a) denotes the transition rates 
between the two states and is constructed from the rates 
of the elementary processes. H will in general be non- 
hermitian. Since H is a stochastic matrix, its columns 
add to zero. The left ground state {tpL\ is hence given by 

a 

with ground state energy Eq = 0, since {iPl\H = 0; the 
right ground state |V'/e) is problem-dependent. All other 
ei genvalues Ej of H have a non-negat i ve rea l part KeEi > 
ijAlcaraz a/.l.ll994llvan Kampe^ Il997|) . 
Since the formal solution of Eq. p62ll is 

|P(t))=e-^*|P(0)) (165) 

the system evolves towards its steady state |P((X))). Let 
r := infi KeEi for i ^ 0. If F > 0, the approach towards 
the steady state is characterized by a finite relaxation 
time T = 1/r, but if F = 0, that approach is algebraic. 
This situation is quite analogous to non-critical phases 



(t ^ 0) and critical points (r = oo), respectively, which 
may arise in equilibrium quantum Hamiltonians. The big 
advantage of DMRG is that the steady-state behavior is 
adressed directly and that all initial configurations are 
inherently averaged over in this approach; the price to 
be paid are restrictions to relatively small system sizes of 
currently L '--^ 100 due to inherent numerical instabilities 
in nonhermitian large sparse eigenvalue solvers. 

Standard DMRG as in one-dimensional T = quan- 
tum problems is now used to calculate both left and right 
lowest eigenstates. Steady state expectation values such 
as density profiles or steady state two-point correlations 
are given by expressions 

(71,) = {^L\nM/{^^L\i^R), (166) 

where the subscripts L, R serve to remind that two dis- 
tinct left and right ground states are employed in the 
calculation. The gap allows to extract the relaxation 
time on a finite lattice and using standard finite-size scal- 
ing theory. The Buli rsch-Stoer tra nsformation (BST) ex- 
trapolation scheme ijHenkel and S chiitz. 1988) has been 
found very useful to extract critical exponents from the 
finite-size s caling functions for de nsity profiles and gaps. 
In this way ICarlon et a/1 l|l999b^ determined both bulk 
and surface critical exponents of the contact process and 
found their values to be compatible with the directed 
percolation class, as expected. 

In cases where the right ground state is also trivially 
known (e.g. the empty state), DMRG calculations can be 
made more stable by shifting this trivially known eigen- 
state pair to some high "energy" in Hilbert space by 
adding a term EsMftlip r){4'l\ to the Hamiltonian, turn- 
ing the non-trivial first excitation into the more easily 
accessible ground state. 

For a correct choice of the density matrix it is crucial 
always to target the ground state, even if it shifted away; 
otherwise it will reappear due to numeric al inaccuracies 
in the representation of the shift operator l|Carlon et all 
l2001bD . Moreover, the trivial left ground state must 
be targeted, although it contains no information, to 
get a good representation of the right ei genstates that 
are joined in a biorthonormality relation l|Carlon et all 
1999b). For nonhermitean Hamiltonians, there is also 
a choice between symmetric and nonsymmetric density 
matrices. It was found empirically ijCarlon et a/.l.ll999hj) 
that the symmetric choice 

p = TrE^i\ibR){^R\ + \^L){^L\) (167) 

was most efficient, as opposed to quantum TMRG, where 
the nonsymmetric choice is to be preferred. This is prob- 
ably due to numerical stability: a nonsymmetric density- 
matrix makes this numerically subtle DMRG variant 
even less sta ble. The same observation was made by 
ISenthil et al\ (|l999) in the context of SUSY chains. 

The method ha s been applied to the ASEP model by 
iNagv et al\ l|2002() . A series of papers has been able 



to produce impressive results on reptating polymers ex- 
posed to the drag of an external field such as in elec- 
tropho resis in the framework of the Rubinstein-Duke 
model ( Barkcma and Carlon', '2003*; 'Carlon et al!, '2001al 
12002: Pacfiens...2003.; Paefiens and Schiitz....2QQ2.) . DMRG 
has been able to show for the renewal time (longest re- 
laxation time) r a crossover from an (experimentally ob- 
served) behavior r ~ ]\[3.3±o.i gj^ort polymers to the 
theoretically predicted r ~ for long polymers. 

Another field of application which is under active 
debate are the properties of the one-dimensional pair- 
contact process with single-particle diffusion (PCPD), 
which is characterized by the reactions A A 00, 
0A A, ^A0 AAA. While e arly field-theoretic stud- 
ies (" Howard and Taiibeil Il997^ showed that its steady- 
state transition is not in the directed percolation 
class, t he first quant i tative data came from a DMRG 
study l)Carlon et all l2001lJ) and suggested again a 
universality class different from directed percolation. 
At present, despite further e xtensive DMRG stud- 
ies jBarkema and CarlonL 120031: iHenkel and Schollwdckl 
I200lir and additional simulational and analytical stud- 
ies, no consensus on the critical behavior of the one- 
dimensional PCPD is reached yet. It may be safely stated 
that DMRG has sparked a major controversial debate in 
this field. While the raw data as such are not disputed, 
it is the extrapolations to the infinite-size limit that are 
highly controversial, such that the answer to this issue is 
well outside the scope of this review. 



B. Stochastic transfer matrices 

iKemoer et al ] l|200lD have made use of the pseudo- 
Hamiltonian formulation of stochastic systems as out- 
lined above to apply the formalism of quantum TMRG. 
The time evolution of Eq. Hlt)5|) and the Boltzmann oper- 
ator e~^^ are formally identical, such that in both cases 
the transfer matrix construction outlined in Sec. IVIII.CI 
can be carried out. In the thermodynamic limit N —f oo 
all physical information is encoded in the left and right 
eigenvectors to the largest eigenvalue, which are approx- 
imately represented using the TMRG machinery. 

In the stochastic TMRG the local transfer matrices are 
direct evolutions in a time interval At obtained from the 
discretized time evolution, 

\Pit + At))^e-^^^P{t)). (168) 

Otherwise, the formalism outlined in Sec. IVIII.CI carries 
over identically, with some minor, but important modifi- 
cations. Boundary conditions in time arc open, since pe- 
riodic boundary conditions would set the future equal to 
the past. Moreover, stochastic TMRG only works with a 
symmetric density-matrix constructed from left and right 
eigenvectors: Enss and Schollwock ( 200^ have shown 
that the open boundary conditions imply that the unsym- 
metric choice of the density matrix, p = Tte\iPr){'iI'l\, 
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FIG. 31 TMRG applied to the causal light cone of a stochas- 
tic Hamiltonian evolving in real time: four corner transfe r 
matrices make up the light cone. From lKemoer 
Reprinted with permission. 

which was superior in quantum TMRG, has no meaning- 
ful information because it has only one non-zero eigen- 
value and that for stochastic TMRG the use of the sym- 
metric choice, p = {1/2)Tie[\iPl){'(Pl\ + is 
mandatory. 

The advantage of the stochastic TMRG is that it works 
in the thermodynamic limit and that, although it can 
only deal with comparatively short time scales of some 
hundred time steps, it automatically averages over all ini- 
tial conditions and is hence bias-free; it can be seen as 
complementary to the method of the last section. Unfor- 
tunately, it has also been shown by Enss and SchoUwocQ 
l)200lj) that there is a critical loss of precision strongly 
limiting the number L of possible time steps; it is a prop- 
erty of probability-conserving stochastic transfer matri- 
ces that the norm of the left and right eigenvectors to 
the largest eigenvalue diverges exponentially in L, while 
biorthonormality (iPlIiPb) = 1 should hold exactly as 
well, which cannot be ensured by finite computer preci- 
sion. 

Bui lding on the observation by lEnss and SchoUwockl 
l)200lj) that the local physics at some place and at time t 
is determined by pas t events in a finite-sized light-cone, 
iKemper et al\ l|2003|) have applied a variant of the cor- 
ner transfer matrix algorithm to the light-cone (Fig. l31|l . 
reporting that the numerical loss of precision problem 
now occurs at times several orders of magnitude larger, 
greatly enhancing the applicability of stochast i c TMR G. 
This method has been applied bv lEnss et all l|2004tl to 
study scaling functions for ageing phenomena in systems 
without detailed balance. 



C. Time-dependent DIVIRG 

So far, all physical properties discussed in this review 
and obtained via DMRG have been either true equilib- 
rium quantities, static or dynamic, or steady-state quan- 
tities. However, time-dependent phenomena in strongly 
correlated systems are more and more coming to the fore- 
front of interest. On the one hand, in technological appli- 
cations such as envisaged in nanoelectronics, it will be of 
great interest to fully understand the time-dependent re- 
sponse of quantum many-body systems to external time- 
dependent perturbations and to calculate transport far 
from equilibrium. On the other hand, the recent mas- 
tery of storing ultracold bosonic atoms in a magnetic 
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trap superimposed by an optical lattice has allowed to 
drive, at will, by time-dependent variations of the opti- 
cal lattice strength, quantum phase transitions from the 
superfluid (metallic) to Mott insulating regime, which is 
one of the key p hase transitions in s trongly correlated 
electron systems l|Greiner et all ^QQ^ . 

The fundamental difficulty can be seen considering the 
time-evolution of a quantum state {ipit = 0)) under the 
action of some time-independent Hamiltonian = 
Enl'ipn)- If the eigenstates \ipn) are known, expanding 
|V'(i = 0)) = X^n^nlV'n) leads to the well-known time 
evolution 



|'0(i)) = ^c„exp(-i£'„t)|-0n) 



(169) 



where the modulus of the expansion coefficients of \ipit)) 
is time-independent. A sensible Hilbert space trunca- 
tion is then given by a projection onto the large-modulus 
eigenstates. In strongly correlated systems, however, we 
usually have no good knowledge of the eigenstates. In- 
stead, one uses some orthonormal basis with unknown 
eigenbasis expansion, |m) = '^j^amn\ipn) ■ The time evo- 
lution of the state \tp{t = 0)) = J2m ^r?x(0)|m) then reads 



representations of operators in the final block bases: 

= [Hes -Eo + Vesimm)- (171) 

The initial condition is obviously to take |'0(O)) — |?Ao) 
obtained by the preliminary DMRG run. Forward inte- 
gration can be carried out by step-size adaptive methods 
such as Runge-Kutta integration based on the infinitesi- 
mal time evolution operator 



m+At)) = ii-iH{t)At)m))- 



(172) 



As an application, Cazalilla and Marston have con- 
sidered a quantum dot weakly coupled to noninteracting 
leads of spinless fermions, where time-dependency is in- 
troduced through a time-dependent chemical potential 
coupling to the number of particles left and right of the 
dot, 



V{t) = -dMt)NR^StiLit)NL 



(173) 



which is switched on smoothly at t = 0. Setting S^l = 
—Sfifj, one may gauge-transform this chemical potential 
away into a time-dependent complex hopping from and 
to the dot, 



|m) 



m), 



tq{t) = tq exp 



(170) 

where the modulus of the expansion coefficients is time- 
dependent. For a general orthonormal basis, Hilbert 
space truncation at one fixed time (i.e. t = 0) will there- 
fore not ensure a reliable approximation of the time evo- 
lution. Also, energy differences matter in time evolution. 
The sometimes justified hope that DMRG yields a good 
approximation to the low-energy Hamiltonian is hence of 
limited use. 

All time-evolution schemes for DMRG so far follow one 
of two different strategies. Static Hilbert space DMRG 
methods try to enlarge the truncated Hilbert space op- 
timally to approximate \4'{t = 0)) such that it is big 
enough to accommodate (i.e. maintain large overlap with 
the exact result) \4'{t)) for a sufficiently long time to 
a very good approximation. More recently, adaptive 
Hilbert space DMRG methods keep the size of the trun- 
cated Hilbert space fixed, but try to change it as time 
evolves such that it also accommodates \i'{t)) to a very 
good approximation. 

Sta tic time- dependent DMRG. ICazalilla and MarstonI 
l)2Q02j) were the first to exploit DMRG to systematically 
calculate quantum many-body effects out of equilibrium. 
After applying a standard DMRG calculation to the 
Hamiltonian H{t = 0), the time-dependent Schrodinger 
equation is numerically integrated forward in time, build- 
ing an effective Hcs{t) = Hcs{0) -\-Vcs{t) , where _ffcff(0) is 
taken as the last superblock Hamiltonian approximating 
H{0). Vcsit) as an approximation to V is built using the 



dt'5nL{t') 



(174) 



The current is then given by evaluating the imaginary 
part of the local hopping expectation value. Obviously, 
in a finite system currents will not stay at some steady 
state value, but go to zero on a time scale of the inverse 
system size, when lead depletion has occurred. 

In this approach the hope is that an effective Hamil- 
tonian obtained by targeting the ground state of the 
t — Q Hamiltonian is capable to catch the states that 
will be visited by the time-dependent Hamiltonian dur- 
ing time evolution. Cazalilla and Marston argue that on 
short time scales the perturbation can only access few 
excited states which are well repre sented in t he effective 
Hamiltonian. With one exception l)Luo et alV 2003) this 
seems borne out in their main application, the transport 
through a quantum dot; in many other applications, this 
approach is not sufficient. 

As on e way out of this dilemma, it has been demon- 
strated ijLuo et alV .2003) that using a density matrix 
that is given by a superposition of states |V'(^i)) at various 



times of the evolution, p = Y^^l[)'^i\'^{'^i)){'^i^i)\ 
^ tti = 1 for the determination of the reduced Hilbert 
space is much more precise, whereas simply increasing M 
is not suffcient. Of course, these states are not known ini- 
tially; it was proposed to start from a small DMRG sys- 
tem, evolve it in time, take these state vectors, use them 
in the density matrix to determine the reduced Hilbert 
space, and then to move on to the next larger DMRG 
system where the procedure is repeated, which is very 
time-consuming. 



with 
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It is important to note that the simple Runge-Kutta 
approach is not even unitary, and should be improved 
bv using e.g. the un itarv Crank-Nicholson time-evolution 
l|Dalev ad 120041) . 

Instead of considering time evolution as a differential 
equation, one may also consider the time-evolution op- 
erator exp(— ii?t), which avoids the numerical delicacies 
of the Schrodinger equation. Schmitteckert ( 2004) com- 
putes the transport through a small interacting nanos- 
tructure using this approach. To this end, he splits the 
problem into two parts: By obtaining a relatively large 
number of low-lying eigenstates exactly (within time- 
independent DMRG precision), one can calculate their 
time evolution exactly. For the subspace orthogonal to 
these eigenstates, he uses an efficient implementation of 
the matrix exponential {ipit + At)) = exp{~iH At)\ip{t)) 
using a Krylov subspace approximation; the reduced 
Hilbert space is determined by finite-system sweeps, con- 
currently targeting all \ip(ti)) at fixed times ti up to the 
time presently reached in the calculation. The limita- 
tion to this algorithm comes from the fact that to tar- 
get more and more states as time evolves, the effective 
Hilbert space has to be chosen increasingly large. 

Adaptive time- dependent DMRG. Time-dependent 
DMRG using adaptive Hilbert spaces has first been pro- 
pos ed in essentially ide ntical form bv lDalev et all l)2004|) 
and lWhite and Feiguia (2Q04J- Both approaches are ef- 
ficient implementations of an algorithm for the classi- 
cal simulation of the time evolution of weakly entangled 
quantum states invented by Ivid al (2003^ 2004) [time- 
evolving block-decimation (TEBD) algorithm]. The 
TEBD algorithm was originally formulated in the ma- 
trix product state language. For simplicity, I will ex- 
plain the algorithm in its DMRG context; a detailed 
discussion of the very strong connection between adap- 
tive time-dependent D MRG and the orig inal simulation 
algorithm is given by 'Dal ev et al\ l|2004 ^ : it turns out 
that DMRG naturally attaches good quantum numbers 
to state spaces used by the TEBD algorithm, allowing 
for the usual drastic increases in performance due to the 
use of good quantum numbers. 

Time evolution in the adaptive time-dependent DMRG 
is generated using a Trotter-Suzuki decomposition as 
discussed for quantum TMRG (Sec. IVIII.C|) . Assum- 
ing nearest-neighbour interactions only, the decompo- 
sition reads H = Hi + H2- Hi = "^^i h2i-i and 
H2 = J2i!Ii hi is the local Hamiltonian linking sites i 
and i + 1; neighboring local Hamiltonians will in general 
not commute, hence [Hi,H2] 7^ 0. However, all terms 
in Hi or H2 commute. The first order Trotter decom- 
position of the infinitesimal time-evolution operator then 
reads 

exp(-iFAi) = exp(-iFiAi) exp(-iff2Ai) + O(Ai^). 

(175) 

Expanding Hi and H2 into the local Hamiltonians, one 
infinitesimal time step t t -\- At may be carried 
out by performing the local time-evolution on (say) all 



even bonds first and then all odd bonds. In DMRG, a 
bond can be time-evolved exactly if it is formed by the 
two explicit sites typically occurring in DMRG states. 
We may therefore carry out an exact time-evolution by 
performing one finite-system sweep forward and back- 
ward through an entire one-dimensional chain, with time- 
evolutions on all even bonds on the forward sweep and 
all odd bonds on the backward sweep, at the price of 
the Trotter error of 0{At'^). This procedure necessitates 
that lipit)) is available in the right block bases, which is 
ensured by carrying out the reduced basis transforma- 
tions on \ip{t)) that in standar d DMRG form the basis 
of White's pr ediction method llWhit3. ll996bD . The de- 
cisive idea of IVidall 1I2OO.I I2004ir was now to carry out 
a new Schmidt decomposition and make a new choice 
of the most relevant block basis states for \ip{t)) after 
each local bond update. Therefore, as the quantum state 
changes in time, so do the block bases such that an opti- 
mal representation of the time-evolving state is ensured. 
Choosing new block bases changes the effective Hilbert 
space, hence the name adaptive time-dependent DMRG. 

An appealing feature of this algorithm is that it can be 
very easily implemented in existing finite-system DMRG. 
One uses standard finite-system DMRG to generate a 
high-precision initial state |'!/'(0)) and continues to run 
finite-system sweeps, one for each infinitesimal time step, 
merely replacing the large sparse-matrix diagonalization 
at each step of the sweep by local bond updates for the 
odd and even bonds, respectively. 

Errors are reduced by using the second order Trotter 
decomposition 



+ 0{At^). (176) 



As At~^ steps are needed to reach a fixed time t, the 
overall error is then of order O(Ai^). One efficient way of 
implementing is to group the half-time steps of two subse- 
quent infinitesimal evolution steps together, i.e. carry out 
the standard first-order procedure, but calculate expec- 
tation values as averages of expectations values measured 
before and after a exp{—iH2At) step. 

Further errors are produced by the reduced basis trans- 
formations and the sequential update of the representa- 
tion of block states in the Schmidt decomposition after 
each bond update during an infinitesimal time step. Ap- 
plications of adaptive time-dependent DMRG so far for 
the t ime-dependent Bose- Hubbard model llDalev et al. , 
2004 ) and spin-1 Heisenberg chains l|White and Feiguiii . 
2004|) have demonstrated that it allows to access large 
time scales very reliably. The errors mentioned can be 
well controlled by increasing M and show significant ac- 
cumulation only after relatively long times. 

Time- dependent correlations. Time-dependent DMRG 
also yields an alternative approach to time-dependent 
correlations which have been evaluated in frequency 
space in Sec. IIVI Alternatively, one may calculate ex- 
pressions as 



\ciit)c,m^) 



(177) 
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constructing both \tp) and \(j>) = Cj\ip) using standard 
DMRG (targeting both), and calculating both \tp{t)) and 
\(j){t)) using time-dependent DMRG. The desired corre- 
lator is then simply given as 



m)\cl\m 



(178) 



and can be calculated for all i and t as time proceeds. 
Frequency-momentum space is reached by Fourier trans- 
formation. Finite system-sizes and edge effects impose 
physical constraints on the largest times and distances 
|i — j'l or minimal frequencies and wave vectors accessi- 
ble, but unlike dynamical DMRG, one run of the time- 
dependent algorithm is sufficient to cover the entire ac- 
cessib le frequency-momentum range. I White and FeiguinI 
l|2004l) report a very successful calculation of the one- 
magnon dispersion relation in a 5 = 1 Heisenberg antifer- 
romagnetic chain using adaptive time-dependent DMRG, 
but similar calculations are in principle feasible in all 
time-dependent methods. 



D. Time-evolution and dissipation at T > 

While our discussion so far has been restricted to 
time evolution at T = 0, it is possible to generalize 
to the (possibly dissipative) time evolution at T > 
l|Verstraete et al l l2004atlZm)lak and Vidail.l2004j) . Both 
proposals are in the spirit of adaptive Hilbert space meth- 
ods. One prepares a T = cx) completely mixed state Poc, 
from which a finite-temperature mixed state exp(— /3_ff) 
at /3~^ — T > is created by imaginary-time evolution. 



prohibitive, but it can be decomposed into a sum of d 
tensor products of A-maps as 

d 

Miioiu'i] = A,[cr,a,] (g) A*[a[a,]. (181) 

ai = l 

In general, d < NsHcM'^, but it will be important to 
realise that for thermal states d — NsHc only. At T = 
0, one recovers the standard matrix p roduct state, i.e. 
d = 1. In order to actually simulate p, IVerstraete et al\ 
l)2004a|) consider the purification 



IV'AfPDo) ^ Tr 



{cr}{a} 



n 



A,\aiai 



(182) 



such that p — Tt:ii\^Jjmpdo){4'mpdo\- Here, ancilla state 
spaces {|ai)} of dimension d have been introduced. 

In this form, a completely mixed state is obtained from 
matrices ^i[criai] cx l-Sa-i.ai, where M may be 1 and nor- 
malization has been ignored. This shows d = A'sito- This 
state is now subjected to infinitesimal evolution in imag- 
inary time, e~^^. As it acts on a only, the dimension of 
the ancilla state spaces need not be increased. Of course, 
for T = the state may be efficiently prepared using 
standard methods. 

The imaginary-time evolution is carried out after a 
Trotter decomposition into infinitesimal time steps on 
bonds. The local bond evolution operator J7i,i+i is con- 
veniently decomposed into a sum of du tensor products 
of on-site evolution operators, 



eM-pH) - (e--^)^poo(e-"^)^ 



(179) 



where (3 = 2Nt and t 0. The infinitesimal-time evo- 
lution operator e~'^^ is Trotter decomposed into locally 
acting bond evolution operators. The finite-temperature 
state is then evolved in real time using Hamiltonian or 
Lindbladian dynamics. The differences reside essentially 
in the representation of mixed states and the truncation 
p rocedure. 

IVerstraete e t al 1 ('20 04a|) consider the T = matrix 
product state of Eq. (|T7|l . where the Ai[ai] are interpreted 
(see Sec. IIII.A|I as maps from the tensor product of two 
auxiliary states (dimension M^) to a physical state space 
of dimension NgUc- This can be generalized to finite- 
temperature matrix product density operators 



{(T}{(T'} 



(180) 



Here, the M[CTiCr^] now are maps from the tensor product 
of four auxiliary state spaces (two left, two right of the 
physical site) of dimension to the A'g^;jj,-dimensional 
local density-operator state space. The most general case 
allowed by quantum mechanics is for M to be a com- 
pletely positive map. The dimension of M seems to be 



Ui,i+i — 

k=l 



(183) 



du is typically small, say 2 to 4. Applying now the 
time evolution at, say, all odd bonds exactly, the aux- 
iliary state spaces are enlarged from dimension M to 
Mdu- One now has to find the optimal approximation 
\'ip{t + At)) to \ip{t + At)) using auxiliary state spaces 
of dimension M only. Hence, the state spaces of di- 
mension Mdu must be truncated optimally to minimize 
IWipit + At)) - Itpit + At))\\. If one uses the matrices 
composing the state at t as initial guess and keeps all 
v4-matrices but one fixed, one obtains a linear equation 
system for this A; sweeping through all ^-matrices sev- 
eral times is sufficient to reach the fixed point which is the 
variational optimum. As temperature is lowered, M will 
be increased to maintain a desired precision. Once the 
thermal state is constructed, real-time evolutions gov- 
erned by Hamiltonians can be calculated similarly. In 
the case of Lindbladian time evolutions, they are carried 
out directly on states of the form of Eq. 1] 180(1 , which is 
formally equivalent to a matrix product state on a local 
N^if^^ dimensional state space. 



The approach of IZwolak and Vidail l)2004(l is based 
on the observation that local (density) operators 
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^Paa'\o'){(j'\ can be represented as A'^gito x iVgitc hermi- 
tian matrices. They now reinterpret the matrix coeffi- 
cients as the coefficients of a local "superket" defined on 
a A'g^;^j,-dimensional local state space. Any globally acting 
(density) operator is now a global superket expanded as a 
linear combination of tensor products of local superkets, 

P= ^ ... ^ c^i...TfJ^i...ai). (184) 

CT1=1 ffl, = l 

Here, \ai) — |(Tj)|cr-), the state space of the local su- 
perket. We have now reached complete formal equiv- 
alence between the description of a pure state and a 
mixed stat e, such that the TEBD algorithm jVidal 
2003V2004") or its DMRG incarnation toalev eA gi.l . l2004i 
White and Feieuin. .2004) can be applied to the time- 
evolution of this state. 

Choosing the local totally mixed state as one basis 
state of the local superket state space, say ct^ = 1, the 
global totally mixed state has only one non-zero expan- 
sion coeffcient, Ci . i = 1 upon suitable normalization. 
This shows that the T = oo state, as in the previous 
algorithm, takes a particularly simple form in this ex- 
pansion; it can also be brought very easily to some ma- 
trix product form corresponding to a DMRG system con- 
sisting of two blocks and two sites. The time-evolution 
operators, whether of Hamiltonian or Lindbladian ori- 
gin, now map from operators to operators and are hence 
referred to as "superoperators" . If one represents op- 
erators as superkets, these superoperators are formally 
equivalent to standard operators. As in the T — case, 
therefore, they are now Trotter decomposed, applied in 
infinitesimal time steps to all bonds sequentially, with 
Schmidt decompositions being carried out to determine 
block (super)state spaces. The number of states kept is 
determined from the acceptable truncation error. Imagi- 
nary time evolutions starting from the totally mixed state 
generate thermal states of a given temperature, that may 
then be subjected to some evolution in real time. It turns 
out again that the number of states to be kept grows with 
inverse temperature. 

For the simulation of dissipative systems both ap- 
proaches are effectively i dentical but fo r the optimal 
approximation of Vcrstrac te et al\ l)2004aj) replacing the 
somewhat less precise truncation of the approach of 
IZwolak and Vidail l)2004j) (see the discussion of the preci- 
sion achievable for various DMRG setups in Sec. lIII.A|) . 

X. OUTLOOK 

What lies ahead for DMRG? On the one hand, DMRG 
is quickly becoming a standard tool routinely applied to 
study all low-energy properties of strongly correlated one- 
dimensional quantum systems which I presume will be 
done using black box DMRG codes. On the other hand, 
there are several axes of DMRG research which I expect 



to be quite fruitful in the future: DMRG might emerge 
as a very powerful tool in quantum chemistry in the near 
future; its potential for time-dependent problems is far 
from being understood. Last but not least I feel that 
two-dimensional model Hamiltonians will increasingly be 
within reach of DMRG methods, at least for moderate 
system sizes that are still beyond the possibilities of exact 
diagonalization and quantum Monte Carlo. 

The strong link of DMRG to quantum information 
theory has only recently been started to be exploited 
- DMRG variants to carry out complex calculations of 
interest in quantum information theory might emerge, 
while DMRG itself could be applied in a more focused 
manner using the new insights about its intrinsic ratio- 
nale. In that sense, DMRG would be at the forefront of a 
growing entanglement between condensed matter physics 
and quantum computation. 
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